{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#本章需导入的模块\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "plt.rcParams['font.sans-serif']=['SimHei']  #解决中文显示乱码问题\n",
    "plt.rcParams['axes.unicode_minus']=False\n",
    "import warnings\n",
    "warnings.filterwarnings(action = 'ignore')\n",
    "from sklearn.metrics import confusion_matrix,f1_score,roc_curve, auc, precision_recall_curve,accuracy_score\n",
    "from sklearn.model_selection import train_test_split,KFold,LeaveOneOut,LeavePOut # 数据集划分方法\n",
    "from sklearn.model_selection import cross_val_score,cross_validate # 计算交叉验证下的测试误差\n",
    "from sklearn import preprocessing\n",
    "import sklearn.linear_model as LM\n",
    "from sklearn import neighbors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0    1204\n",
      "1     892\n",
      "Name: 有无污染, dtype: int64\n",
      "截距项:-4.858429\n",
      "回归系数: [[0.05260358 0.01852681]]\n",
      "优势比[[1.05401173 1.0186995 ]]\n",
      "预测结果： [0 1 0 ... 0 0 0]\n",
      "总的错判率：0.153149\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAECCAYAAAARlssoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO29fXRcZ3no+3v3nk9p9GU5kR0nRnZCnaSHJBA1QBJCapKSQHqAXrhwKfesttyEtjSc0kUpXT2UU0oL7aLnclZW6Vq+JZweFqe5CW3IKfQ0IRgSEpKmNiU63GCf2LGdSLFkS9ZImu89e7/3jz3v1sxoRhpJM9LM6PmtJVuzZ/aed4+k53mfb6W1RhAEQRAArK1egCAIgtA+iFIQBEEQAkQpCIIgCAGiFARBEIQAUQqCIAhCQGirF9AIO3fu1KOjo1u9DEEQhI7i6NGjM1rri9ZyTkcohdHRUY4cObLVyxAEQegolFJn1nqOuI8EQRCEAFEKgiAIQoAoBUEQBCGgI2IKtXAch4mJCXK53FYvZcuIxWJceumlhMPhrV6KIAhdQscqhYmJCfr6+hgdHUUptdXL2XS01szOzjIxMcG+ffu2ejmCIHQJHes+yuVyDA8Pb0uFAKCUYnh4eFtbSoIgNJ+OVQrAtlUIhu1+/0L7cPz4cY4fP77VyxCaQEcrBUEQBKG5dGxMYTuQTqfp7e3d6mUIQl2MdZBKpSoeHzhwYMvWJGwMsRQ2yOHDh3nLW97Cbbfdxi233MJtt90WfP/8888D8Mwzz5DNZtFa4zhOcK7Wmnw+Hzx+9NFH+c3f/E3OnTsHwN13383Jkyc394YEQdjWbBtL4dS0ww+PFZhZdNnZZ3PjlRH2jWw8lfPgwYMcPHgQgBtvvJEnn3xy2WsmJib4zGc+w3333cev//qvE41GAV8p9Pb28vd///dkMhls2+ZDH/oQsViMgwcPMj4+zszMDBMTE7zwwgsbXqsgNBtjEYiF0D1sC6VwatrhkeeyxCOKoV6LVM7jkeeyvOsGmqIYDKHQ8o/T8zxuvfVWbrnlFnbs2MH3vve9iue11riui+M4fPe732VwcJBrrrmGb37zm/ziL/4ijz32GHfddVfT1igIgrAS20Ip/PBYgXhE0RP1vWU9UQV4/PBYYUNKYX5+nne/+93Bzv8nP/kJBw8eJBwOo5QilUrx+c9/nt/93d/ly1/+Mu973/tIp9OcPXuWq6++GvCVxh133MEnP/lJPvShD/Hxj3+cn/u5n+NLX/oSn/70pwHfX/v1r3+dY8eO8cd//Mcb+zAEoQWIhdA9bAulMLPoMtRbGT6JRxQzi+6GrjswMBDs/F999VU+/vGPc9FFF/GJT3yC8lbfDzzwALFYjCeffJIHHniA48eP85nPfKbiWplMhj/90z/lK1/5Cn/xF3/BHXfcwcDAAACf/OQnefnll7n77rs3tF5BEITV2BZKYWefTSrnlSwEn2xBs7PPbtp7PPzww7ztbW9jfHx82XPxeJzPfvaz/OVf/iXnzp3j8ccfD5TJ5OQkL774Ik888QTnz5/nN37jN9Ba8+ijj6K1Zv/+/ViWheM47N+/n7179zZtzYIgCNVsC6Vw45URHnkuC3jEI4psQZMtaH7hukhTrv/qq6/y13/91zzxxBMVSmFiYoLdu3fzD//wD/T39wPwsY99jI997GMAJJNJ3v72twNw5513cuedd+K6Lh/5yEfIZrP09PQwODjIV7/6VelvJAjCprAtUlL3jYR51w1xEjGLubRHImbxrhviTQkynz17lrvuuovPf/7z9Pf3Y1kW58+fB+Azn/kMjz32GF/72td4z3ves+zcV155hT179gB+wPmpp57i9ttv57rrruNXf/VXee9738sb3/hG3vzmN3Po0CFeeeWVDa9XEARhJbaFpQC+YmhmphHAyy+/zNve9ja+8IUvcMcddwDw7ne/m0984hO4rsuOHTu48cYb2bt3L2NjYxXn3n333Xzve98LAsff+ta3OHToEPfddx8/+7M/y/3330+xWOTee+/lne98J4cOHWJqaorLLrusqfcgCIJQjtJab/UaVmVsbExXj+P86U9/ylVXXbVFK1piYWEhcA2thUwmQzwe33D/onb5HARBaD+UUke11mOrv3KJbWMptIr1KASAnp6eJq9EEARh42yLmIIgCMJW0WkdZEUpCIIgtCFbpUzEfdREZmZm2LlzZ9OuJ11SBaFz6dQOsmIpNJF3vOMdvPzyy8uOS5dUQRAaxVgIqVSKVCq16RbDtrIUWqGpb7rpJvr6+sjn8xw/fpx77rkneK5QKHD33XcTCoWkS6ogbDM6tYNsU5WCUioEvFT6ArgXeC/wDuA5rfVHS6/7o+pjncrTTz8NwEc+8hHuvvtuPvjBD1Y873kes7Oz0iVVEISG2Gpl0mxL4Rrgb7XWvweglLoeuBm4AfhDpdRtwFz1Ma31401eRwWt9u09/PDDPPzww5w8eZL7778f8Duo3njjjXz84x/nAx/4gHRJFYRtSqdYCIZmK4U3AXcppX4e+J/AceDvtNZaKfUocCcwX+PYMqWglLoHuAdo6yZwjzzyCJ/61KcYHR3lt3/7t4PjL7zwAtPT04yOjkqXVEEQ1sxWKZNmK4V/AW7TWp9VSv1XII6vGAAuACNAEThZdWwZWutDwCHwK5o3sqhWmWNHjhzh0KFDPPTQQ9xzzz3EYrHguUgkElQrS5dUQRA6hWYrhXGttUmnOQKE8RUDQAI/2ylV41hHMjY2xre//W1Onz7NyZMn+cIXvhA8Nz8/z1vf+lYA6ZIqCELH0GyB/DWl1LVKKRt4N9CLHz8AuBY4DRytcWxTOHDgQMtMsje+8Y08/vjjwdcXv/jFwFKQLqmCIHQKzVYKnwW+BvwYeAb4HPB6pdR/Bj4F/C3wVI1jHU2xWKx4/Nxzz/HhD3+Ym266iUKhULdL6nve8x7e9773AX6X1D/7sz/jvvvu47d+67dwHIdCocC9997Lgw8+yEsvvcTU1NSm3ZMgCNuTlndJVUrFgXcCP9Jav1Tv2Eq0c5fU9SJdUgVBaDVt2SVVa50FvrHasXVee8NCdatoRpfUTmh7LghCZ9GxQd5YLMbs7Oy2FYxaa2ZnZysyngRBEDZKx7a5uPTSS5mYmAhGX25HYrEYl1566VYvQxCELqJjlUI4HGbfvn1bvQxBEISuomPdR4IgCELzEaUgCIIgBIhSEARBEAJEKQiCIAgBohQEQRC2gK2awbwaohQEQRCEgI5NSRUEQTB00sjLVg/92ihiKQiCIAgBYikIgtCxtPuuuxZbPYN5NcRSEARBEALEUhAEoWNp9133SrTrWsVSEARBEALEUhAEoeNp1113JyKWgiAIghAgSkEQuoitrJJt1wpdYW2IUhAEQRACJKYgCF3AVubrd2KtgFAfsRQEQRCEALEUBKELWG++fjN29a2qFVjP9cRK2ThiKQiCIAgBSmu91WtYlbGxMX3kyJGtXoYgdA3VcYBEIgG0xw57PWtr5/vZSpRSR7XWY2s5RywFoSOQdEdB2BwkpiAI25B27hm0nrW18/10GqIUhLZG0h0FYXNpiVJQSo0A/6S1fr1S6ivA1cC3tdafKz2/7JggCJtPOyvX9aytne+nU2iVpfBFIK6U+iXA1lq/WSl1v1LqtcDrqo9prV9s0TqEDkfcAoKwuTQ90KyUOgikgSngVuDB0lOPATfXOVbrOvcopY4opY6cP3++2csUBEEQatBUpaCUigCfBj5VOtQLTJa+vwCM1Dm2DK31Ia31mNZ67KKLLmrmMoUO5MCBA2IldCj1Mscko6w9abal8Cngy1rrZOlxCoiXvk+U3q/WMWENyB+TIAitotkxhduAg0qpjwLXAXuBV4BngWuB48AEvsuo/JggCF3G8ePHmZycJJvNsmvXrmUbGckoa0+aqhS01reY75VS3wf+LfADpdQlwJ3AmwBd45jQAJKe2X3Iz1BoN1pWp6C1vhVAKXUrcDvw51rr+XrHBEHoHoyyGxgYIJfLMT8/TyKRqFB+ohDbk5YXr2mt51jKNqp7TFgdSc/sHjrF6mvXdQmtQyqaBUFoOuUbmGoLofo1QnshSqEDkT+mzqfdrb5OsWSE5iNKQRCEliFKpPMQpSAIW0i7Cs12t2SE1iGFY4IgCEKAWAqCsIU0Yyfeyt28WAjbD7EUBKFLkXYownoQS0EQtoBmZPdsZoaQxBa2D6IUBKHLkHRSYSOIUhCELaAZ2T31rtFMl5EomO2HKAVB6DIknVTYCKIUBGELaYbAbqXQFwWz/RClIAhdighwYT2IUhCajuwqm0crPsv1XFN+ltsHqVMQBEEQAsRSEJqGZKo0j1Z8lvLzERpBLAVBEAQhQCwFoWlIpsrqNPrZNPpZHj9+nMnJSfbs2dO0awrbG7EUBEEQhACltd7qNazK2NiYPnLkyFYvQxDWTbU/P5FIAOvfrRsLIZfLkcvliEajxOPxhiwGYfuglDqqtR5byzliKQhCGyAdTYV2QWIKgrAJNNuff+DAAQ4cOLCmmIIgNIIoBUHYIBsR9JImKrQbohQEYRNptrA3FoMgNAtRCoKwTpqxy1+rW0ksCaHVSKBZEARBCBBLQRDWSTODx41aCBJ7EFpNS5SCUmoHcD3wr1rrmVa8hyC0K/UmoZUL8Or001rPrSTwJycnV32NIKyHpruPlFJDwLeAG4DvKaUuUkp9RSn1jFLqP5S9btkxQWgVK9UBbLRGoDzYa1JEG2VycrKh15v3SCQSJBIJ9uzZw549e9a9ZkGoRysshWuA39FaP1tSEAcBW2v9ZqXU/Uqp1wKvqz6mtX6xBWsRhE2jvMo4lUpx+PBhAAYGBiqe37NnDydPngRgbm5u2XMruYiMAim/ZvVrBGEjrKoUlFK9wBgwgm9ZnAb+Wdfpj6G1fqJ03i341sIO4MHS048BNwOvr3GsQikope4B7gHYu3fvGm5p81jrH6T8AW8+K/nim+mnN0I9m82Sz+eZnp4mnU4zPDxc8brZ2VlgSRk4jgP4wn52dnbV3b953qxZEJrNikpBKfUrwPuAJ4A5IAHcAfxnpdRtWuvFOucp4P2lczRg7OMLwBuA3hrHKtBaHwIOgd/7aC03JQjNphGFYXb5U1NTxGIxLr/88grlU+vcXC4XnGuqkld6vXQ6FVrNapbCh7XWb6k+qJT6EvB24Bu1TipZER9VSv0x8F7g/yk9lcC3NlJAvOpYx7DWHWY7ZY5sN2GykhA13xs3z/XXX1/zGuWunUbexzSmK39f4/Y5ceIEAPG4/+sfi8WWrUsQtpLVlMKcUuoP8V09k/iC/C3AbcAXap2glPo94KzW+r8Cg6XX3Qw8C1wLHAcmahwTythuwrtdqY4TNGoxrLbjb0TB1FqLeV5+L4RWsZpS+CXgXuDP8F0+i8BR4B1a66k65xwCHlRK/V/AT4BvAk8qpS4B7gTehO9S+kHVsY5hrSZ8O5j87WStNJtG7qXWc+Y848KplUpaHSeYn59fVaCbbKbyz9pYGr5ndeV1CcJWsqJS0FoXgf+79NUQWus54PbyY0qpW0vH/lxrPV/vWKeyXgG7Uv56NwrvTqQ6TrAV3Ujld0LYTDalormkKB5c7VinUW0BlLOSD3sraLa10kxFuNE1bFRY1hs0VStO0Mi1a33W1eeJQBfaldWyj54EeoAFQOG7fRR+LPlg65fX3qxXKK10Xju4mrqZ9X6uKymEVv+s5HdC2ExWsxTeB/wX4P1a64XWL6fzMVkm7WrqN8tCaKYiXC8bFZaN7t7Xu8aVUkpbTbv93gmdw2oxhWml1P8BeJu0no5ivUKpkfPkj7m5rKSU1ttHaLN9/fI7IWwGq7mP7i/7PvgW3330ay1cV8fS7aZ+KxXhWmhmemZ1TcFam9O1ExKUFjbKau6jHvxeRg7wPH466o+AMy1eV0fRTPeC0BpqKaXq1NH5+fmK167nmgYRxkKnspr76AMQtMJ+P/Ar+Ompx4GrWr24TqbbhcFWKcJW7ISnpvySG2MNN9qcztBObay73VIVWs9q7qMH8CuOC/iWwt8Cn8RviicIHUctoVndeXStlBezNdKATwS10M6s5j4aAv4V30o4ga8grsVPTZWYgrDpNHMnbJRBrQrlRlpV1FIAq/VI2ixE8QjrZTWlYAPfBp4BbtBa/2rrlyQIraNckGez2UCIHzhwIGiMtxGhXkuZ1Gp7YZ7batppLUJ7sJpSCGmtvw6glPp+65cjCI2xESE2Pj5OoVAgHo+jlOLpp59mfHw8mH3QqOtnpYp2QehUVlMKFyulPoifhmq+B0Br/d9aurJtiuzcWouJIczOzpLL5YLW1YuLi/T29pJMJjl79iyDg4NcccUVG3qf6u/b6WcrqatCPVZTCv8v8Noa38vQmy2gkYHw25G1VlUPDAyQy+U4duwYc3Nz9Pb2AjAxMUEul2P37t3E4/E1u37MiE3Ddv+5CJ3Jaimpf7RZC9nutMPObaub5W3We0xPTwMQCoVwXReATCZDLpfD87xASdRzC603BbWdlIRsLIR6bEqXVGFjVCuMWgPhoXv+sBu9n/HxcSYnJ2t+DtUCvbzy2Qj1m266CYDnn38erTWJRIJisbjsvEbXVK/bqiB0EqIU2oSt3Lk1y0rZDGun/D0cx2F2dpZsNsuuXbsaOu/AgQM1s4vi8TiDg4MsLCyQTqeJxWLL1l2v2WEn0y0bCaF5iFLoAKoVhpkl3O4WwlrdLI0qFTMeMxaLEY1GyWQyzM/Pc/311wcxAOPfN7t3MxO5uj355OQko6Oj7NmzhxMnThAOhxkeHl6mOFa6B5mVIHQTohTajK3w5zfLSqm+TiKRWNd1alFdgZxOp+np6WHXrl1Bz6JaJJNJoLKFhblOdYbQ5OQk8XicgweXRoXIrAthuyFKoYPolB1pdduItc6yNjETYxFVU76LTyQSK1pOZi1mDrOhnlWykcK1dv15CMJaEKXQJTTDn98soWYEq1nLRii/r2PHjjE+Pk5vby8jIyPMz88zOTm54rrNWqpdWbXiAeXVyCs1xBPhL3QzohSEprNeN0t5HYF5bCqNazWrq56DsFIlcjabrbmmZriCxJ0kdBOiFLqE9Qi51V7bDsLOtJ6Yn59naGgosBBMYLl8naYFdq1K5NWyk8x1TBB7ZGQkeN+V5jMLQrchSkFoGeUVwY0I1eqsoEQiwcDAAFNTU2QyGRzHCYR1rfNqZTutpiybYSFIqwihmxClsAHaUQisxUKo525Zy4CZzWDXrl0kEolAUVQHoM36TDDZBKrLs4hWovyeBwYGyGazzM/Pk0gkGr6GIHQLohRaSLu0e+i0grj1pn+Wu3xqXbPZSIqq0I2IUlgHne42qCXMyhu/maBuq/3pa51KVu+56lTWRlNhy3+O5fd8xRVXdMzPspxO+z0U2hNRCi2gFUqjXtbMSu+xlnVsdCRlNevZRdeLCawlLrEVtOu6BGE9iFJYByvttDdTQExOTjbU96ce5Wut1dahloJpxv2VK6vp6engPmBJMZSvo9axWhj//1qL5Tq9BqHTLVehvWi6UlBKDQAP4I/yTOPPd/4r4Grg21rrz5Ve95XqY91CM5XGan/wJvhafd3y4qxabqDq4i1zfeOPX0uMYiUhVO/eTeonwOzsbMUktGolUK0wRNgJQutohaXwy8B/0lp/Ryn1V8AHAFtr/Wal1P1KqdcCr6s+prV+sQVraSnVvnjY2C6t0XPLu3XmcrmmCstqi6GWUqoeTl+vWrg6VbSadDrN8PAwl19+eaAkTAO7XC7H3Nwcs7OzzMzMEI/HCYfDNa+/UTpdyUjAW2gmTVcKWusvlz28CPgQ8KXS48eAm4HXAw9WHatQCkqpe4B7APbu3dvsZbYc8wc6OTlZt+d/I6z2Bz8yMlLReM687sSJE8BSI7jy52BJyJtzzf8rrWtqaop4PE4ulyOVSgWB3epqYYN53lgEZgzm8PAw2WyWQqHA7Owse/bsYc+ePYHyMOmwWmvi8Thzc3NEIpGmxTsEQahPy2IKSqk3A0PAacBsFS8AbwB6axyrQGt9CDgEMDY21rbTS+oJbbOjNkJwJdbqE96KnWEikSAej1fUL8zOzgJg2zbgD70B2L9/PwCnT58GwLIs5ubm8DyPcDjMmTNngqE2tm0HloexPowFVp0FVK44at2z+TzqjcXc6OfV7jvxdl2X0Fm0RCkopXYA9wH/G/A7QLz0VAKwgFSNY11B+U49mUyilKro0Fmv8+dqrNWyWO14edVweUrn5ORk3YKt8oZxANdccw2wJIRjsRjgj7tMJpMUi0V6e3tJp9O4rkskEmF4eDiIITiOg23bFcqg3vsKgrA5tCLQHAEeAn5fa31GKXUU3z30LHAtcByYqHGso6m2EJLJJI7jAH5P/1wuR19fX0Pn1gocr3aeeVyr++d6qHb9GL9/uYvHvP/4+DiLi4vcdNNNHDhwgMOHD5PL5YJZx47joLUmn88Hn0U4HKZQKKx6b+X3ZeIOqVRqxV373NxcxeONZudIdo+wnWiFpfBhfHfQHyil/gD4KvB/KqUuAe4E3gRo4AdVx7qCckFhfO3xeDxwtWxW2moj2UMDAwOkUikeeughYMkNdPjwYU6fPs3o6GjN883O/fDhw+zZs6fCPXb8+HFmZ2cpFosUi0UcxyEej2NZViBUXdclHA4TCoXwPE+yigShjWhFoPmv8FNQA5RS/x24HfhzrfV86dit1cc6mVq57iamEI/H6e3tZW5urq4ArBbUa2lfkUqlmJqaqujwud5CtYWFBZLJZJDZNDc3RzqdJp1OB2mjU1NT5HI5stksp0+fZteuXdi2zYkTJ5ieniabzaKUCoLcfX19RKNR+vv7GRsb4+mnn6avr4/e3t5gMlo9F5GxEGBJyZq+RPU+m8HBwYrHG40pSHaPsJ3YlOI1rfUcS9lGdY+1E2ttqlaNERyzs7OcPn2aXC6H67qk0+nAathohXN5amgymSSTyTA3N8fc3BxjY2N1z60ePGNiA+Pj43iex+joaNCGGnwX0uLiItFolIWFBRYWFgDwPA/XdZmYmCASiRAOh8lms2SzWUKhUHAM/PRT407q6+sLYgu5XK5tW1OLEhC2I1LRvE7qFX/VshjGx8cpFouAn8VTnY3UaA+g6ufKA78miGtY6fzq7CGjJBYXFwGCTqHpdJp8Ps/u3bu55pprmJycDDKKQqEQruti23bgBopEIlx88cWcP3+eQqFAJpOhv78/ONes16x5dnaWaDRa0zKqlVo7Pz/Prl27KobrVMceylnt8VoR5SBsB0QpVGEsBOPWaNRiqDXgxQg24zoqFotorTeUTVNeuGbWNzs7y9mzZ5mbm8N13VXXbdZqAr1Hjx4FIBqNAn720MmTJwNFViwWSaVSgRIZHh5GKUUqlaK3t5fdu3cHysJkG4VCIRzHYWhoqG7mU29vL7Ztc+zYMSKRSGCxbDUSWBa2M6IU1ki1wDBFX/G4n2FbT3CUF5jVu1Y94VP+OlPBfPbsWQYHB8lkMoC/czdCHJasgOprHDhwIKgnML738owiIIgZxGKx4LpHjx5laGiIm266KTjH1GAcPHiwIjtoYWGBXC6HZVkUCoUV5zVMT0/jui7Dw8M1P7vqKWqNFAOK8BaE9SNKoQqzs15pp22EMhC4bUzqZiPN1apdJGvBVDCbALYRmidPnmRiYgLHcbj44otXLJgzzxnFYa6xZ88exsfHyWazuK6L4zg4joPnefT09FTcf3mxWfk9mZRa27YJh8MMDQ0tC64bN9Lhw4dJp9P09PQwMDCw5TvytfR0EoRuRZTCGikPIIMvYMurfNdzrUaFT/n0sd7e3oqiOLOrN7UR5VPIZmdn6e3tRWtdkaH005/+FCAQ+E8//TTJZJKRkRHy+Tz5fB6tdRAwdhynZsdSs/4jR44EytIU7Z04cYJEIkFfX19Fj6jytNfyLKhq91i1hVU+9wDWXwwoCEJtRCnUoZaFYFwfU1NTzM/Po7WuKMCqFytYy06znoKo1VQumUwGbhzw20kopTh37hyO47B3716mpqbIZrMMDg6STqcZGhpibm4uqLYGX8CmUilCoVAQVygXxpdddhmFQgGtdaBsjAuqep2RSATwC9ZMDKVYLFIoFCoUEsDo6Gjgelrpc6p2bzUbiSEIwhKiFNZJPB4nk8mQTqfrvqbRyV+rPZ9KpTh79mwQjDUzBsLhcIWgNMVnJlU0m83S19dHsVhkdnYW27Yr3Er5fB7wrQrHcejt7eXChQuk02ksy8K2bWzbZnFxMYhXpNNpEokEWmvS6XTgZtuzZw+Dg4MMDg4Sj8eDwHO5JZBIJBgfH2d+fr4iJlCeWrua9WTagDdiIYhwF4S1I0phDZSnU9YKtMLKLabrxRCqg6/lx2BJaIO/Q3/66acpFotEIhFOnDgRNKozbqL9+/dz9OhRkskk4XCYdDpNJpPBsqxAuQBBm2rzHufPnw/aV5vZBgZTY2DaVAwODuI4zrK01nL6+vqWxR1qWTz16hSqXUnz8/NBV9VmIjEEQVhClMI6MAqhUCiQzWaXKYDqFtPVyqC6oriesDNdQY01Eo/HOXv2LAsLCwwPD9Pf3193jYlEInh/s8s3u/vp6Wmi0WigAIx7JxqNBpXIJqYwOjoaKBETe8jn88zMzBCJROjp6Qkyr0zAeqWBP41ORys/Pj09HazfrGGldiHiDhKE9SNKYR2YzJvy4HJ5IVl1i+labRqAICg7ODgYCLvqKWnm3EKhEAybcV2X3t5eXNclHo8HgrZ8x2tGdP70pz9FKUUoFApqC4rFYuAC8jwPx3ECZWC+PM8D4Ny5c8H7DA4OcuWVVzI/P8/p06cJh8Ps2rWrZrptOfXqFFbb8Vd/Bq2epyBKQxBEKayJU9MOPzxWYGZxNzv7bC6LnWFk0F4mTKpbTBuhaapzTX2AsQBWmrFsFNDjjz9e4c83wWRThQy+sDUBYIN5DxNADofDwffGgojFYiilsCwLz/OwbZtYLIbruhQKBWzb5qqrrgKWlJZJX52fn6/r319ph76W1ha15kevhLiDBGH9iFJokFPTDo88lyUeUQz1WqRyHv/8cp4rdoc4MpVmZtFlZ2S1mbEAAB8RSURBVJ/NjVfuZ9+I3+9nOulybKJINjrqK5GdLiOD9rJrm53w9ddfv6yN9tmzZxkaGsJ1XSzLHztRLBbp6+sjlUqRzWaDbqXgt6ro6+ujUCiQTqeDHb/pTLpv3z4mJiaCim2TOmpZFrt37wb8QrhsNltRDHf27Fny+XzQ/ru8O2qjwne9bp2N1HUIgrA2RCk0yA+PFYhHFD1RXzD3RBWpvv08M+lx+S4vUBSPPJflXTf45/zo3KXEE4qhiAqUyBt/Jhoohmr3yEoYARyJRDh16hShUIiRkRHS6TSnT5/mxz/+MaFQiFwuFzSaKw8UO46DZVkkk0kuXLiA53nB80ZxGPfW7t27eeWVVwiHw4HFYFpX9Pb2MjU1xezsLD09PezatavujOZq4b9R1tvdVBCExhGl0CAziy620rx4wSVX0MQiioKjcbWuUBTg8cNjfu2C62omL3jB68ORUb77kmJHn8XOPpvp5AuMDNoMDAwws+By6O/GWci62Ja/gy/kIoSLCeL9kYqGdydOnGBmZgbwBbrWOphiVj6nQCmFbdvLsowikQiJRCKIaYBf47CwsEChUODcuXNBgVwkEglcVNlslmQyGdQsmHoHQz3hb5RGo0FmQRC2DlEKq2DiCGfOF8kVoCcCsYjCKWqSGc1AT+Xr4xHFzKJLJudxIa0JWaDQzCxoHBdAk5k9yfkwnI5YvOXqKNmMy/OnC0TjESwFU3Medn6a/lienLPAT19K0x9aZG5ujiuuuCJID02n0xSLxcCtBL5wHxgYCIrMTAGZCR6bwHEotPSjV0rhum5gMSwuLgZWgcG4jcLhMAMDA+zatYuXXnoJoKIFBiy3GFYLRAuC0D6IUqiBUQQvn3eYS2sGehQFBxwXFrLgaU00YmEpjad9F0x65iSZgseM9xpcD3KORgF2WJHKg+tBLH8KgIIFRQ8Wcx6P/jjHzr4w0XiEHbtey4tnHUL2SZQNGTdKfygK0UFc3Qt6mtOnT+O6LkNDQywuLuJ5HpZlBV+ZTAbHcYJGdkbQG2vBdV1c1w0a6VmWhdY6+AKCFNXp6Wmuu+46AGZmZoJpaufOnSMejwczEZLJJENDQ8s+x+qsIbEQBKH9EaVQRXlAOVuAoquZmPWFZSQERRdSOXBcD1vBYlYzfqYAi45vCSQ0V+yyOTHlknegUNTYFngaIo6fb+/aCTQQclPk8nAqO0IiCqF+Fyd5Eis/TbGYw/MAvUDc9cgXHaxQjkQiQaFQqBiCY9pll1chx2KxoGDNKAYgSDkNh8MVgWRYUhy2bQdprKdPnyabzQbWSG9vb1AJbYLb5W4og1EAjfSFWo+yEAUjCK1BlEIV5QHlnFOk6IJlgeeBpfzviy4UimBbEMmewsmD7aWxFVjZUyzOWLgli0EDduYUMUBp36UTzZ8DwLP94jOtYTEHL/1/TxEqJvF0gZC7iLb78HSRVKZAxDnLPEVisRhaa1zXDVJLjdBXShGJ+PGHfD6/TOgDFW6icoVh23bw+kKhQCQSob+/P0hpNcrGZCiZSmbjqorFYkClkDZKo7r2olnIbGdBaD6iFKqYWXQZ6vV3xbGIIpXThGxfKbierxDAF+RFF8LKVxZuaTPuaDg751GM+QrBYKyEYmgntucrBye0s/LNNRTtQSydQ3sZbF3A04CbAc/BQ3P+/PkguJzNZtFaBz2PbNsOAtKvvPJK8LpaGIshEokEE9Q8z8PzvCDeUCgUgnGa5R1ZV6I68+ill15icXGxZqFa9WsbGWhUa7YEiGIQhGYhSqGKnX02qZxHT1QxMmAxu+j5wt/24wAG2/KFfi6yDw+IuaeWntS+krAtiORPBcrB9nLkrV4ysdoCzNI5LHcB28sCfqwi5Po7cg+wYEWXj4kVOI4TZAiZmEGw7pICMTMTgCDt1ASk+/r6gn5L+/btI5fLkc1myefzQSqqbdtB+qvp/wRLQtu0pjDtM8xzzRDeptuqadMhikEQmocohSpuvDLCI89lAY++uMWuQYtXL3hYFijtKweAeBjSBV9Y10LjWxcJdxqvkARdQJfiCoXwCLnoviDwDP5xXyFksL00WoVRXoH67+Dv9o3Aj0QiFZZBeY1C+essywpcT5ZlBbOWyykUCoRCIS6++OLAXRWNRoPg9EqUz6WGpYl0tQS3+d5YCI0EpMtrO7TWDbXZEAShcUQpVLFvJMy7bqDUzsJl/0iYn/83FmfOe4yfKeB6moIDWafkNsLPKurJHadoD+LaCcLFmeCxVzyH5WWxvAzaiqBCF6PcaXLsC94zF93nu5dUBMgAGsvLAQ4KI9yXu4HKBT34isHMSTDPmayinp4ewuFw0Ouov7+fK664omJs59mzZ8nlcuzevZvdu3cHDfkymUzQshsqexfV252bmgrjdmpmZ9PqKW9iIQhC8xClUIN9I+GgVYXhZuDGad+KmJx1KeZ1EEcAsD3fL+7aS7vWkOsXdmkVBgVKF3GKUAiPEMufCuIMAOHiOSwvg+0uonDw3UcaGlAKJlU0FAoFSsBYCsZdNDw8TDKZZHBwEK114Aoy6aWmDsGct1JL7EaoHmu6kuBeT1GbKAJBaA2iFNaAsSLu+8dFih7EC6cI5af9wDAWtpukJ5fCUxFy0f3YboqoM4GrIigVAXwBbZSByUaKONNYXgbX6sF2FwGFZ/VgeSl0STGo2ksC/JhAoVCoqGI29Qv79vkWSTgcDobu5PP5oM7A9Eq65ppr6OnpIZlMsnv37kAppFIpBgYGmJycZHJykoMHD644f7qaVk1LA1EMgtAKRCmskULyJeL5FFm1j5DlWwNKF1AUUdrFdhfQyqZo96NVDNfyLQfX6sMtpaCG3CSeFUPpAp7yZxW49iBFezCIKfgKJAR4aBVFaT/OUY1Sip6ensBqKC9W8zyPubk54vE4H/zgBwF/537u3Dk8z2PPnj1BvyNYmiYHBDOgzfMmNlCLRmIAjSBCXhC2HlEKwFMvZPnO83kWsh79cYvbr41y89XxslbZpgNqaf6wC66CBXsfsZi/0w8Xz6HtGKrk5tEqRiE8wkLiTfSnniXkJtHKz+XPxA5ge2kSmaOAzXziprLV+DaB5WUAP9OonkIAX/jPzMwEFgL4isEMwDEYwZ3L5QKrYXx8PDinfNZzebM708rCtPdup6KxdlqLIHQL214pPPVClm88kyUcUvRGFZmC5hvPZDk/73JmxiMeUeQLHj8af5Hn/hUSoTSuA1H8zKFc1HfPxPMnSlf0x2ZGC6cJF88F71O0/RkKhbA/tN52U4CNZ8WCa/RljpaUh4XGKjmNKjODalGuEEw8oVgskk6niUQiKKV4/PHHg0E5O3furEhHhaW5DeBnDmUyGRKJRMVMZ6jMIiqvGZienl41+CwIQvvTEqWglBoBvqG1fotSKgz8PbAD+IrW+v5ax1qxjkb4zvN+y4a8o8nk/doCS8H3Xyhw9aVhFjIep855RD2/VUU6XzrRqryOv7PXaGUDFrZOYXspIs40oeI5PLsfDUQLZ4gVJlC6gJ9lpNkx/0/YXgrbXfCDzDoEFAOrYzVM3yMzQc1kI5kCtZ6eHjzPo6+vL6gZMKMzy9t2T05Oks1mSafTFAoFxsfHg6yjdmpuJ+M2BaF1NF0pKKWGgL8BTIvNe4GjWuv/qJT6R6XUQ8Dd1ce01ov1rtlKLqQ8PA+UBWjIO361sgYcV3PmvIvrQTrs7+bLawvKKYZ2YHkZlM6jlY1n+fGDQngES+ew3QXCxfN+/YF2gvOULhAuzvhxiZKiUBRLAebVMbUGoVAI27Yr0lONAN+7d2+Q4fPQQw8BtQWomTtdLBaDucyTk5McOHCgbnM78NtYaK0rnhMBLQidSSssBRd4P/BI6fGtwKdK3z8JjNU59r3yiyil7gHuAV+otQrL8quPlfZjBeWi+PnTy3sHVRPEC7DQ2Fg6j9JFcpHXACbeMFMKOCvAQiv/Y7c80zpCVfxvFIIOzqiNaVoXiUQYHBxkdHSU2dnZIKPIdDFdiWrhPTk5ydTUFIODg1x++eXBceNaWq253WYg4zYFoXU0XSlorRegoqK2FzA+igvASJ1j1dc5BBwCGBsba2zbvA6Gei2m5jyKulYlwJJlYPz+hmjhDOC3rog4r+LrQt91BBYhdwGAcPEllM5iWX7gWGkXpcuVTcnVo/OomsFkU6+wnHA4jFKKwcFB+vr6gp2+SVHNZrNcc801wJIANcHnlQSqmbhWq/q41nnl7iUR0ILQ2WxGoDkFxIF5IFF6XOvYplCdUTQyaBO24fT5SoFsW0tN7kxdQbViiDhTuPYgrtVbChArtIqhtEPEmUSrKJbO48cOMvhBaKtM+NuAMU88loR/yZdFBE+F0DpN+WRn4yoyPYhGR0eDnfzJkycBljXCM+4f0/a6Xr+gVtYVNBtRQILQfDZDKRzFLwj+BnAt8GydYy3l1LTD//hRlv91tkhPRHHJDr/x3VzKxbIUYRuUAjt9Cso6n9peGqVzRJzpigpk213wBb6XwtJ+kBkslM6i8HzXj7ZYSiU11kGo5B7SpUCyLr0mhMYrva7kisIF7WCriP8aXQwa1pnGeIODgxX3GYvFgvYUsLSLN43ryrONDNWB25WoJ4hFQDeGWFRCu7MZSuFvgH9USr0FuBr4Z3zXUfWxlmEG55yb94iHFRrFy+ddXnNxiOE+G9fTLPRaXFj0fGeOBldDT+64X2BmxXyR7y21j7a9DEZQW16BpbiAF8QCNIXg9aYm2VOmfUZ55ECjVQTLyxOkoCoLDxulXVy7B9tLYykraGsdjUYJh8PE4/GKFtL1GB4eXqYsqjGdTY2VIQJMELYfLVMKWutbS/+fUUrdjm8Z/KHW2gVqHWsZZnBO0dVEQr5FAIrppMsVu0LMpTWXRk6RynmlamJwLT95yhfWOVw7UsoOWkJpD4s8WllYOge6ULMdRalv6VI4WRfRKlQqSvNjEErn8KxeoAdwUdoJ1mB5RWx8RRCLxUgkEoyOjgIE7aPN98lkMig+Az/byHEcYrEY2WyWTCazrPbAYKwIk0UkNA9JoxU6hU0pXtNavwo8uNqxVmEG58QiCqfoD80JWf4c5WxB4xQ1J19diimEijPYKoVn+UIyXEyhvRzFoE2FH0RWulgWNF4eEFYVz+iSbeDhWTE8FcX2TPWzH2x2QjuX6hW0i+Vl8aw4TmgnsVCYwQGboaGhYH6B6VYKS2Mvk8nkip+FsRhqUX28EYHVTsKtndYiCJ3KtqhoNoNzRgYszpx3Ma6bkOXXKZyZdnEjS3UItkpRCI8EVkO4OIPtpYKqZEvPgApRtAdLnVB9l5FWPaAzNa0FU4jmKRvLy6KUbwn4gWg/jJyJHfCL3awktpvE9lJYukjeHkTZDpCvEOrVwm9+fp6BgQF27drF/Pw8AO985zvrtrqu3r02WpgmwnftSBqt0ClsC6VgBufEI4q9Oy1endOksh6JuMUrM+6yRhJmCI6pQXCtRKAQQm4yiC0YSwI8X+jXcR8toVDa9V9f5jHzrBiWlyZWOEPYmcIqFbIpXQga6kX6RnjrLfuWCZPyQTXlbiNDuTXRKI0ILHPNdnCHiGtGEJrHtlAK5YNzXp7xShPKIJX1cKo0QnnaadSZKHU9Vf6sAy/lVyNrD0URy5unvFGdWrFPke9EsrQTuJF84V8kH7406JNk6QKWl8XPTlIoXOLFCYrZBLCv7tXLexcBXH/99asOw6neva5GufDNZrPBWMyRkWVlJkIdRFEJ7c62UApAMDTn/HMuFxYVntZ1d/Wx/Cn6Mkex3Xk0FpaXw/JypfoCha5IJ62MJJRmoWFSTksTlCFobmcT1CuUAtchN0nIXcAJXey30rZ0aRRnEVSI3kSC8I6rarogqnfJxm0ES72MUqlU03fPZgxmO0w/E9eMIDSPbaMUnnohy0PPZMk7pfoDCyJ27df6Tewu+D2M7F6sop8qqkuZQiYmocoKzvwjpU6lhFD4LiYVKBGFxg6G53hWGE/5rh5/jsJSjYCnoliWQnkeKjJAz8VXs3P3axu+V9PaIpfLkc/nmZqaIh6Pr1qYtpr7RYSvIHQ/20IpPPVClgeeylIo6y7hepAtLH9tf+pZos4EfksKh1CxgO8iMrt9cxEbj3BgMWhlo5UfG9BWBE0EpZ0gduCV+h05dj9RL11KZ3XwrBiunSBjX4lr9fqT2qwEni5gWREIDVJwNDdeGVnRd17dxdRkI01PTwcFba0Q4u2kGNppLYLQqWwLpfCtI9llze5qEcufoif3QtCawq8wLk0yI1SyDDz8VhUaz4ri2IN+CwsVwrUSRAtplHYo2juwdLZUkAaF8CUU7UG/HsFN4lk9lQHsUhosugglt5KnevAiI1yxO8S+kTDHV8g2Le9iauoM5ufnG1IIa7UAVnteLAlB6Fy2hVJIZvxZCKvhVzCbquKlmIMGULbfUxsbz+pFaQfP6vFrF+x+IsUZAPLh3dheBqUdNLY/glO7eFaMUCnNFMBTETKxJaFp0mBd29/p226KaBj2XfEzZMJ+HUN5phH4wWRDvS6mnTz0RpSLIGw+20IpeLUnWVbQn3qWaOHlijgBUIojmC6mvsWA9tAqjGslgrGbZtSm0jm/3XVpt58PX4rSS+0x0EV0KfYQcaaDSWwmDdZ0ZXWtXqywP9NhZrF2VlO50Ky129+oBbBWoSypoYLQ+XS9Ujg17aBUaZNfh1j+VKkIrfJFS4+WGttpFSpN5LGDYrOIMx0M0tEq4g/oCe3EKykKozCKdgzXihEuXqh4n+ruq+aYa8PxySJXX+ZnThkha9xD1X2MylNQm8VKvZJahSgXQdg6ul4p/N0zmVVdRz05X+hoFa4qQLNKuUMWRXuIkJsM2loUQkPBqwrhEWKl+QrgN8tTgBeK+T2RoEJBaBXCtfsDK6EcoyDClj8AyHFhLuVyanppWtvU1BSwNLOivM6g3F20XuFafl55s71GLQ0R4oLQuXS1UnjqhSwvnl29117RHiTqTPjZQhXPhNClbCNTiQwKrcLko/srXqmBot1PPvKaQMmUt8owze0Mxl1ksKqsGdNpdbAHhvtsfniswC+/1Rey1aMxVxqVuV5MYZppuLeZFoMoF0HYOrpaKXzn+XxDAeaFxJvYMf9PhDhf0dHUn5Vso7H8amZ8wa+xiRYmyMQOVOz2bXfBn7mgC3il5nm2m1qmAMpR+BZBUNqmwLb9x/GIYkfun8m9CrmLbgzOqXYPmTTU6iBzo8K1+nnz/+TkJFrroFCtUUSIC0Ln0tVKYS69coQ5lj8VDM4JF6fxexgBJZeRmaesVRjl5fziMxWvaKFdLuxN4NhkENXDnGP6qiogZIOloej6Wam2Da+5OERxSuG4mt19S5V29YR9M3fW5YphqyqWRbkIwubT1UqhkayjcjQRFH6Ngo8qfdkUbd8146ehXsJC4k3Lzq/OIMpF99W1ECh7l6Lnp8wmYmBFFOm8Zkf2WYpTFrn0PK6G+MIzHD5scfDgweD8RoXmahZCvZiDCGVB2H50tVKwrdrHjdCOloLDfnaQjbaiaK+IHzsIl+YeRHDtwWUDdqqpJfzLlUMja80WoCeqCVnQF1dkC5qwrdg9YDHYW+dmyuj2imVBEFpPVyuFVfpYV+C3r84Q9DZSRggrivYg4eI5XLufXOQ1q17LKAGjFFZbom357qOiC3kHrro0xG/eeSewVKhWbiE0i424nSQILAjdSdcqhVPTDtn8yq8xOiMfHqEQHiHiTBN1JrC8NEW7MuXU/N/Irt8oA5N5tJLFoPEVwkCPheNq8o7mzjfEV32PTkOUiCB0Bl2rFL7xw/SqvY5q4VoJPBXxs4xKtQWNKIKV6ItDKrv0uHpwZ28U8kXfbfQzl4SDNt/QPAthJaG8HgtBCssEoTvpWqVw5nz9KHO1kC9/XAiPBLGGlc5ZiXL3kaWgGN/PRUNwNqn94HepJkFr6InC1ZdFyBb8edF3vD628sU7DFEigtBZdK1ScNdhJqykLFbDUv6Xp5fSTI1FMDJgEY8qQiHNYtbDcf0ActhWjAxazKU9dvbZ/MJ1kQoroRk0WyhLYZkgdDddqxQSUUUys7Jm2KhbCPxiM8uCUKktRa5A0GspH9tHNAQXDdjMpT0+cHMPPzxWYGbRZWefzY1XNl8JtBuiRAShs+hapdAX1yQz6zt3LcpCA/09CkuBUwQd0ThFCIcgElLEo35q6c4+m30j4U1XAq0SyiLcBaE76VqlkM5Xh3Objyr9c9nOEFprTk4V2TtsMzVfimdozUCPTbag+YXrIi1dS7vT7UpELCGhW+hapdBMYmG/fgB8RWBC2CHbtxL64xaZvMfPXBKmP25RKDrkHIhFLHYNtoebSISVIAiN0LVKYfTiEBdSzuovXIHhBNw1FmfPcIj/cjjF+QXf8gjht7SO2HDpsE0m75EtaN51Q3zLhb+wuUh2ldBtdK1SuOP1MX700tqUgsIPEgMkYoov/Du/gO3rT6S5ZEeIiwdhOumSczSepwmHLFwPhnqtlmQOCYIgbDZbqhSUUl8Brga+rbX+XDOvvW8kTDwM2RX0QsgC11uKPOjSP5EQXL5rqSvpzKLLUK+FUr6rCEBrzVza49/f1d/MZQsdhmRXCd3G6l3WWoRS6pcAW2v9ZmC/Uuq1zX6Pqy4L0xPxhbxVsgAUfsHYG/aHed1rQlyywyJs+6ml0RD098AlO6yKVhM7+/xgcTkmo0gQBKGb2EpL4VbgwdL3jwE3Ay+aJ5VS9wD3AOzdu3ddb3DH62Mk015FwVhf3OIDN/dUuHpOTTsr1g/ceGWER57LAh7xiAqqj7d7RpGwhFgIQrewlUqhF5gsfX8BeEP5k1rrQ8AhgLGxsXXllu4bCTdUMLZa/cC+kTDvuoGK60gMQRCEbmQrlUIKMD6aBC1yZTWrYGwrCs8EQRA2my2LKQBH8V1GANcCp7duKYIgCAJsraXwTeAHSqlLgDuB5fMtBUEQhE1lyywFrfUCfrD5WeDntdbzW7UWQRAEwWdL6xS01nMsZSAJgiAIW8xWxhQEQRCENkNp3dpOos1AKXUeWD4ObXV2AjNNXk4nIfcv9y/3v33ZCfRqrS9ay0kdoRTWi1LqiNZ6bKvXsVXI/cv9y/3L/a/1PHEfCYIgCAGiFARBEISAblcKh7Z6AVuM3P/2Ru5/e7Ou++/qmIIgCIKwNrrdUhAEQRDWgCgFQehwlFI7lFK3K6V2bvVahM6na5WCUuorSqlnlFL/YavX0mqUUiNKqR+Uvg8rpf5BKfW0UurX6h3rFpRSA0qp/6GUekwp9bBSKlLrZ9+tvw9KqSHgW8ANwPeUUhdtp/s3lP4G/rX0/ba5f6VUSCn1slLq+6Wv1yml/kgp9S9Kqb8se92yY/XoSqWwGVPd2oWSUPgb/PkUAPcCR7XWNwHvVUr11TnWLfwy8J+01r8ATAEfoOpn3+W/D9cAv6O1/hPgUeAg2+v+DV8E4rXutcvv/xrgb7XWt2qtbwUi+N2nbwDOKaVuU0pdX31spQt2pVKg9lS3bsUF3g8slB7fytK9PwmM1TnWFWitv6y1/k7p4UXAh1j+s7+1xrGuQGv9hNb6WaXULfh/9G9nG90/gFLqIJDG3xTcyva6/zcBdymlnivNvH8b8HfazyB6FHgL8NYax+rSrUqheqrbyBaupaVorReqOszWuveu/zyUUm8GhoBX2Gb3r5RS+BuDOUCzje5fKRUBPg18qnRou/3+/wtwm9b6BiCMP7hsQ/ffrUphU6a6tSm17r2rPw+l1A7gPuDX2Ib3r30+CowDN7K97v9TwJe11snS4+328x/XWp8tfX+EJtx/N3045WznqW617r1rP4/STvEh4Pe11mfYfvf/e0qpf1d6OAh8gW10/8BtwEeVUt8HrgN+ke11/19TSl2rlLKBd+NbBRu6/64sXlNK9QM/AL5Laapbtw/xUUp9X2t9q1LqNcA/Ao/j7xrfBFxafUxr7W7ZYpuIUuo3gD8Fni8d+irwO5T97PFdKl35+1BKNHgQiAI/AX4fP260Le6/nJJi+LdU3StdfP9KqX8D/DdAAf8d35X2A3yr4Y7S15nqY1rrU3Wv2Y1KAYI/ltuBJ7XWU1u9ns2kNOL0ZuBR88tf61i3Uutnv51+H+T+t/39x4F3Aj/SWr9U71jd87tVKQiCIAhrp1tjCoIgCMI6EKUgCIIgBIhSEARBEAJCW70AQWg3lFL/Eb8Y7Dzg4FeNJ7XW/7tS6gEgB/x74AHAxq+mfb/WulDjWiHgpdIXwL1a6//Z8psQhHUiloIg1OZPtNa34Ke4zuLnd1P2f3XPpTvqXKeiN40oBKHdEaUgCCszhG81FJRSw/iWQ62eS+fqnF/Rm6ZkOQhC2yJKQRBq8wdKqSfxhfoj+MVx72epSA5Y6rmktX62znWqe9O8o3VLFoSNI0pBEGrzJ1rrW7TWvwzMAz8CfqX0P7Cs51I9qnvTdFPbZqELEaUgCI3xI+DnSv/X6rlUj+reNM+v8FpB2HJEKQhCY5wG/hd+HxmADwNvwHczfV8p9X6l1NVKqc9VnfdZ4GvAj4FntNaPb9aCBWE9SJsLQRAEIUAsBUEQBCFAlIIgCIIQIEpBEARBCBClIAiCIASIUhAEQRACRCkIgiAIAf8/Ghutulw9LNoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "data=pd.read_excel('北京市空气质量数据.xlsx')\n",
    "data=data.replace(0,np.NaN)\n",
    "data=data.dropna()\n",
    "data['有无污染']=data['质量等级'].map({'优':0,'良':0,'轻度污染':1,'中度污染':1,'重度污染':1,'严重污染':1})\n",
    "print(data['有无污染'].value_counts())\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111)\n",
    "flag=(data['有无污染']==0)\n",
    "ax.scatter(data.loc[flag,'PM2.5'],data.loc[flag,'PM10'],c='cornflowerblue',marker='o',label='无污染',alpha=0.6)\n",
    "flag=data['有无污染']==1\n",
    "ax.scatter(data.loc[flag,'PM2.5'],data.loc[flag,'PM10'],c='grey',marker='+',label='有污染',alpha=0.4)\n",
    "ax.set_xlabel('PM2.5')\n",
    "ax.set_ylabel('PM10')\n",
    "plt.legend()\n",
    "\n",
    "##Logistic回归\n",
    "X=data[['PM2.5','PM10']]\n",
    "y=data['有无污染']\n",
    "modelLR=LM.LogisticRegression()\n",
    "modelLR.fit(X,y)\n",
    "print(\"截距项:%f\"%modelLR.intercept_)\n",
    "print(\"回归系数:\",modelLR.coef_)\n",
    "print(\"优势比{0}\".format(np.exp(modelLR.coef_)))\n",
    "yhat=modelLR.predict(X)\n",
    "print(\"预测结果：\",yhat)\n",
    "print(\"总的错判率：%f\"%(1-modelLR.score(X,y)))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "说明：\n",
    "1、这里以空气质量监测数据为例，对是否有污染（二分类输出变量）进行预测。首先对数据进行预处理，将质量等级是优和良的合并为0类（无污染）共计1204天。其余合并为1类（有污染）共计892天。这里只考虑PM2.5和PM10对有无污染的影响，作为输入变量，只有0和1两个取值的有无污染作为输出变量。建立Logistic回归模型。\n",
    "2、首先，modelLR=LM.LogisticRegression()定义modelLR对象为Logistic回归模型；然后，modelLR.fit(X,y)表示基于给出的X和y估计模型参数。其中，X为输入变量（矩阵形式），y为输出变量。\n",
    "3、模型参数的估计值存储在intercept_和.coef_属性中，依次为截距项和回归系数。从回归系数估计值看，PM2.5（系数为0.05）对是否有污染的作用比PM10（系数为0.02）更大。\n",
    "4、modelLR.predict(X)表示将X带入回归方程计算y的预测值。\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "训练误差: 0.0782442748091603\n",
      "混淆矩阵:\n",
      " [[1128   76]\n",
      " [  88  804]]\n",
      "F1-score: 0.90744920993228\n",
      "AUC: 0.982303010890455\n",
      "总正确率 0.9217557251908397\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmMAAAETCAYAAAB6AgEhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd5gUVdbA4d+ZQE5DVoKAoGAA1CGDDBIkSgYBlYy66uoa1rhrWBUD+ClmsqJkFBRQUBElJxEVRJEoCILkMHnO90f1BMYZZhi6pzqc93n6qdBV1acVitO37rlXVBVjjDHGGOOOMLcDMMYYY4wJZZaMGWOMMca4yJIxY4wxxhgXWTJmjDHGGOMiS8aMMcYYY1xkyZgJGCJSzO0YjDEmL8RR1O04jH+yZMz4NRHZLiJXiUg4cFJECmZzXFMRaZxhu42IdPWs259zY4wrRCTJ80PyUuC3cxzXXkSuyLDdS0Ra2v0rNNj/ZOM1IvKUiJwSkYMisk9EHvDsHy4if4jIARG5K8Px14jIZs97z2Vz2TjghKomAymqGp/Ncf8DymXYvhKIEZGhwOsX/u2MMcFORAaJSKyI/Om5Lz12jmOzvN9lIVZVTwHxnldW1woDRgOFMuxuBtQHnhaRh/L2jUygsGTMeNsbqloe50bysIjUA14EYoCmwP9EpI6IRACzgceAakBrEWmXxfUUSMmw/jcicj3wB7BORHp5dp/GSeQmA2s9xxhjTE7mqWoFoC4wUESanuPYrO53mSV7lucaYf0W4CMgTkRae/al3sNeAP4UkYbn8yVMYIlwOwATnFR1l4isBtoDi1T1VwAR+RzoBKzH+cU4z7N/LtAaWCwiG3H+bB4HqgMfiUgCEC4iyz0f0QwoChQDXgJ6A/8GCotIMnAj0BBo5LnOeuBbn39xY0xQUNW/RGQ+0AJYmcOxqfe7y4FNInIU55FkPFDMc98qAFT0rBcCqqpqeRGpBdwJdATGAz+ISBmcH6+VgD7AEaAIsNYHX9X4AUvGjE+ISFUgGigOrMnw1h6clrA44JcM+yfhJFcAicC/VHWpiPwE9FDVvSKSpKrNPa1qiZ5XI6AqMNdzfivgEc+1z6jqQE88V4lIFVX93Sdf2BgTjIT0lvnsD0q/3/3HsysR6O1J0o557luVgeWe9ZrAN55j2wMXAcuAP4GZwECce9gyVX3S8xjzOhEpo6qHvfkFjX+wx5TG2+4WkYPANuBlYAdn95NIAAoDpYBTqTtV9U9V3eHZzPHml3oaMB/n1+O3wD+AQ8BUnEcDLUTkc88v1mc9xxljTI48ydNNwPWe/q4ZX897Djvrfqeqmzz7c3MPS31s+ZaqVgM2A3cBu3HuawA9RORLnJa5h4AyF/zFjF+yljHjbW8AzwC/AwuAGpzdKbUgcAbnl2NaZaSItASqqOoHnl2vi8i5HlMCoKoqIrfgJFrbcVrdhgPrcH5s9AE6AydVdbV3v6oxJgh1FZEDOH22XlHVN7M6SESe4u/3u4xmici5HlOmUhF5BNiH88P1SuBmYAtOy/9NwFBgQ2p3DxN8rGXMeJ2qngEm4rRU7cBJkFJdAuzE6U9RI8P+FjgdZlPdo6rNPcf28Kwne5YxqQeJSCTprV4vALVVdZiqvgvswrmJPQuc9N43NMYEsXmqWlFVL80uEcso0/0uo96e+9Upz7IHcMCzfnOG4yoA9wLX4lSFh6vq3ar6Os5ThDs9r9MX+sWM/7KWMeMrbwDfATcA/xWRy3Faw9oBTwL7gXc8lUNrcDrgP3K+H6KqiUA1z/hjdYHrMrz9K04iVp9zjO9jjDEX6A3gOxF5XFXPK2lS1f3ARZ4BYesDlwGrPf3EdgP/wknUDnk5ZuNHLBkzPqGqu0XkW5xm9oeApTgtsY+r6i8AItIJGAuUB8ar6mee04VcPqYUkTrAIpyWr2XAChGpBozyfN5anMrKFiJSVVWf8tV3NsaEpgz3uwE49zQhl48pRaQN8B5wAFgOfCUiVwNjcH5QVgYaA9eIyO+qOjEfv5rJJ6J6rqFPjMl/IvI9cJ+qLs3m/dRqyqKqesaTZO3xvHcdTt+N/6rqWBGJ8mxHAzGqes4SdWOMuVAicgyor6q7snm/JrDSM7RFBFBeVf/wvNcVJxG7T1U/9vy4/By4GLhGVbfnw1cw+cySMRN0RKSAqiZk2A7DudkdcDEsY4zJlSzuYQWAEqr6l4thGR+yZMwYY4wxxkVWTWmMMcYY4yJLxowxxhhjXOSTakoRqQDMVtUW2bwfiTMpamlgQk7VIWXLltVq1ap5PU5jjP/asGHDX6pazu04vMHuYcaElvO9f3k9GfNUr71H+jyDWbkHZzThp0RkoYjMUtVsB+WsVq0a69ev93aowa1TJ1i40O0ojDl/nn6sIrLb5Ui8xu5hxoSW871/+eIxZTLQFzhxjmNicCZDBWdOwejMB4jICBFZLyLrDx0K4rHuOnUCEe+/LBEzASKBcE5RwO0wjDHGNV5vGVPVEwAicq7DiuLMwwVwBGc6iMzXGYszeB7R0dGBXfLpVitVx46wIPN0acb4j6NHY+nZcyaFCkXwySf9iIiwbqzGmNDj1p3vFM6cWwDFXIzDNzK3duWUiHXs6Dya8fbLEjHjx3bsOErTphP5+utdbNx4gN27j7kSh4hUEJFl53g/UkQ+FZEVIjIku33GGJNXbk2HtAFoDswG6gGrXYrD+7JrBbNWKmPSrFr1O127TufQoTNcdVV5FizoT9WqJfM9jrz2cQWGZ953rn6vAIwcyeFjCZQu7PxGMyZHlSvD4MEQFlztFebvfJ6MicgNwBWq+kaG3e8BC0WkBXAFzkTRgSurBMySL2OyNGPGTwwcOJf4+GTatbuUWbN6U6JEQbfCSe3jOu8cx8SQPol9ah/XrPZ9fa4POvTc/9Hg9ABu5Dde5zMKkHwBYZuQccUV0KSJ21EYH/NZMqaqMZ7lEmBJpvd2i0hbnNax/6pqYN6VrBXMmPPy+ee/cfPNcwC4/fbreOONjq72E7uAPq459nsVkRHACICqVauysdc9/PlhMmOTovm56nXM6a2UO1d7nAltU6bAzp1w8twNriY4uPWYEs+kqDNzPNCfnKsjviVgxuSodevqtG9fkzZtqnP//U1ySoL8RWof1+M4fVxPZbPvLJmLkNpN/g/f3rWPbt1msGzPSRrMLsW8eTdTr17F/PoeJpCsWuUkY926QUQEPP44PPxw7s9XhaQkSEj4+ysqCkqX9l3s5ry5lowFFEvCjMmzo0djAYiKKkxkZDgLFvQnLCwgkrBUWfVxzVO/1wYNKrFu3XC6d5/B2rX7aNZsIu+/350ePer4KHQTsBo2hC++gFjn7w+PPAKLFmWdXKW+EhPP3j6X5cvhkkuc4woVgosv9v13MtkKiInCo6Oj1bUBEzMnYpZ8GZNrO3YcpVOnqVSsWIxFi26hQIHwXJ8rIhtU9W9jEPqCiCxV1Zis+riKyCXAQuBLoCnQGKiced+5ultkvofFxSUxfPinfPDBD4jAjz/eyZVXlvfJdzMB7OhR2LMH6tfP2/kREVCgwNmvvXuzPvbVV+Gmm6BaNasw8YLzvX9Zy1h2LAkz5oJkrJiMiAjjyJFYKlYs5nZYWcpDH9cL6vdaqFAE77/fjbp1y3P8eLwlYiZrUVHO64cfYN++9IQqMvLvSVbmV2Rk1lWYY8c6rWxhYVCwIPzxh7P/vvuc1+DBMPGcMxQaH7CWscysMtKYC5a5YnLmzF6ULFnovK6Rny1jvpbbe9imTQeIiirsyjAfJkQtWwYDBzr901I9+yx07w6nTzuPMK+6ylrLztP53r9s8JLMMreG2eCpxuSaqvL888u4+eY5xMcnc/vt17FgQf/zTsRC0YEDp+jUaSoNGoxjxYo9bodjQkWLFrBjB2ScdvCJJ+DKK51+a3Xrwr//DcmBOehBoLBkLFXqqPmpLAkz5rzNnLmZxx9fggiMGtWWt9/uZFMc5VKBAuHUqVOOgwdP06rVe0ycuNHtkEwoKVsWvvsOataEKlWgdu3090aNgtdecy+2EGB3Sci6f5gx5rz16nUF/fpdxZw5fXjggaaBMnSFXyhdujCffTaAf/6zIYmJKQwd+gn/+tfnJCWluB2aCRXXXAPbtjlFAz//7CRnqb7+GubPh/h49+ILYpaMZUzE7LGkMedtx46jHDx4GoDw8DCmTu1J9+42VENeRESE8dprHRg3rguRkWG8+uoaOnWamjY8iDH56ppr0lvE5s+HLl3g6qvhxRedpM14TWgnY5kTMUvCjDkvq1b9TuPG4+nWbTpxcUluhxM0hg27lq++uo1y5YqwePF2Vq783e2QTKjq1g169Ejf3rbNqcZ86CH3YgpCoZmMpfYPs0TMmDybMeMnWrV6j0OHzlC8eEESE62Drze1aHEJ69YNZ9y4LnTqdJnb4ZhQVbUqzJnjDIFx991O6xjAmjXQpw+0aQPXXQfVqzvDcAwbBtu3w08/QYo9Ys+t0BzaImM/FkvEjDkvqsrIkct5/HFnOC5fzTEZikNb5GTVqt9ZtWov//pXY+uPZ9yxYQNE5/Kv5SWXwIMPwqlT0KuXUxwQImzQ15x06pS+HgCJqDH+JCEhmTvumM+kSd8jAi+/3DaQ5pgMaKdOJdCz50z27z/FDz/8yTvvdKZQodC7hRuXXXstTJ8Ohw8781uWLu20iBUv7rSaHTwIJ044x+7eDffc46w/+ij897/w9NPuxe7HQutvcuY+YsaY8/LBBz8wadL3FC4cwYcf9rCO+vmoWLECjBnTgYED5/Lee5v45ZfDfPxxX7+d1cAEKRHo2zfr91I79Z844Yzmf+aMk5Ct9kzd+swzcNFFcMcd+RNrAAmdPmPWWd+YCzZoUH3uvrsB33wzyBIxF/TqdQUrVgyhatWSrF69lwYNxrFhwx9uh2XM2UqUcKZUmj4dVq48e+ioqVPdi8uPhUYyZomYMXm2evVe/vjjJABhYcLrr3ekQYNKLkcVuurXr8i6dcNp1qwKe/eeoHnzScyatdntsIzJmgh06ADTpjnby5bBe++5G5MfCo1kzBIxY/Jk5szNxMRMpnPnqZw+neB2OMajfPmifPXVbQwZUp/4+CSKFi3gdkjGnFuLFunrgwY5Q2PMnOlaOP4muJOxzFMcWSJmTK44FZPL6Nt3NvHxyTRsWImCBUOri6m/K1gwgvHjb2L9+hF07FgrbX9KihUmGT9UqRJ8/HH69qhRcMstEBfnVFuGuOBNxmyKI2PyJCEhmWHDPuGxx2yOSX8nIlx77UVp20uX7qJBg3Hs3HnUxaiMyUbXrvD8807n/rAwSEyEokWhZMmQbyUL3rurTXFkzHk7ejSWDh0+ZOJEp2LS5pgMHKrKf/7zNd99t58GDcaxdOkut0My5mwizhAX//d/UK+esy8lxXnddZeTnIWo4EzGMo4lZkmYMbk2ffpPLFmykwoVilrFZIAREebP70fHjrU4fDiWtm2n8M47Xhws2xhvWrPGmZD85Zed7b/+OvtpVogJzmTMxhIzJk/uuCOap5+OYc2aYVYxGYBKlizEJ5/czEMPNSUpKYU771zAXXctsKmqjP+JjIQqVWDgwPR9e/e6F4/LgjMZS2WtYsbk6KOPfmb37mOA07ry3/+25JJLSrkclcmr8PAwXnqpLe+/342CBcN566313HzzHLfDMiZr5crBbbc563ff7QwK+8IL8OWX7saVz4IvGcv4iNIYky1V5YUXltOz50w6dZrKmTOh218jGN16az2++WYQlSuX4M47g2KKTxOsmjZNXz9wwOlX1rat03IWHg7z5kFsrHvx5YPgSsZsuiNjciUxMZnhwz/l0Ue/QgQGD65P4cI2dEWwadSoMtu23UObNjXS9qW2ghrjN26/HU6fhmHDoE+f9P179zqd+7t1gyJF0qdVCkLBlYzZ4K7G5OjYsTjat/+QCRM2WsVkCMg4mfjixdupWfN1nn9+Gao2HpnxI0WKwLhxMGOGk4R9/bUzsXhGjz3mTmz5ILiSsVSWiBmTpZ07j9K06QSrmAxRW7YcIjk5hccfX8KAAR8RG2uPpo0fqlQJYmLg6aedQWFTJyb/+muYPNnNyHwmeJIx6ytmTI4+//w3fv75L666qrxVTIag++5rzNy5N1OsWAGmTfuJFi0msXfvCbfDMiZ7BQvCW2+lbw8eDIcPuxePjwRHMmZ9xYzJlTvuiObttzuxfPlgq5gMUTfddDmrVw+lRo0oNmxwBohdvTp0hxQwAaB0aVi5Mn37RPD9gAiOZMz6ihmTJVXl//5vFdu2Ob8kRYQ77oimZMlCLkfmH0RkgoisEpEnsnm/uogsEJFlIjLasy9CRPaIyFLP6+r8jfrCXXlledauHUarVtU4cOAUvXvPIj4+ye2wjMlekyZQtaqzvm6du7H4QHAkY6ksETMmTWrF5P33L6ZLl2kkJNjAnxmJSA8gXFWbADVEpFYWh70I/E9VWwCVRSQGqAtMU9UYz+vH/Ivae8qUKcKiRbdwzz0NmTq1h00Eb/zfnj3O8r//hdq1g6pDf3AlY8YYwKmY7NAhvWJy5MjWFCgQ7nZY/iYGSJ2deDHQPItjLgO+86wfBEoCjYHOIrLW07KWZRYjIiNEZL2IrD906JB3I/eSyMhwxozpQIsWl6Ttmz17C8ePx7kYlTHZePhhZ/nLL85r5Eg4c8bdmLzEkjFjgkxqxeRXX1nFZA6KAvs860eAClkcMxt4UkS6AO2Br4B1QBtVbQhEAll2VFXVsaoararR5cqV83rwvjB//q/06TOLxo0npD3aNsZvPPggvPsuvP9++r7q1d2Lx4ssGTMmiKxevZdGjcbz889/ceWV5axi8txOAYU968XI4n6oqs8CnwHDgPdU9RTwg6ru9xyyHsjq8WZAuvLKclx5ZXm2bv2Lhg3H88UX290OyZh0ZcvCiBFw661wxRXOvoMHYexYd+PygsBPxmxIC2PS/Pjjnxw6dIZ27S5lxYohVjF5bhtIfzRZD9iVzXHfA1WBVzzbU0SknoiEA92ATb4MMj9Vrx7FypVD6Nr18rRH3WPGrLEBYo3/2bgxfX3UKPfi8BKfJGO5qFCKEpGFnv4U717Qh9mQFsakGT78OubM6cP8+f2sYjJnc4FbReQVoA+wWUSezeK4h4BXVDW1c8ozwBScJG2VqgbVjMbFixfko4/68vjjLUhOVu6993OGD//Uqi2NfylQAN5801nfti3gp0ryejKWywqlW4EPVTUaKC4ieZvFNmOrmFVSmhCUmJjMvfd+xpYt6R3Ee/SoQ2SkddbPiaqewOnEvxpopaqbVPVvPyBV9UlVnZJh+ydVrauqV6vq4/kXcf4JCxOeffYGpk/vSeHCESxcuI0jR4J7omYTgPr3T19v0gTq1AE/LZbJiS9axmLIuULpMHCViJQCqgC/Zz4gV5VI1ipmQljqHJNjxqyld+9ZJCenuB1SwFHVo6o6U1UPuB2LP+rb9yqWLRvMvHk3c9FFxd0Ox5izlSoFL72Uvr11K9x7r3vxXABfJGO5qVBaDlwC/BP42XPcWc6rEslaxUyIyTzH5KRJXQkPD/wuoMb/XHfdxWcVgYwcuYzZs7e4GJExGTz0EOzfD8WKOdvTpjmPLQOML+7eOVYoAU8Cd6jqM8BWYLAP4jAmKGVVMdmwoVVMGt9bu3Yfjz22hN69Z/HUU0tJSbGO/cYPVKwIy5alb7do4V4seeSLZCw3FUpRwNWeaqRGgP2NNiYXZs/eQqtW73Ho0Bnatq1hFZMmXzVocDGjRrUlLEx4+ulv6N17FqdOJbgdljFQvz506OCs//mn02IWQHyRjOWmQmkkMBY4DpQGpp33p9iQFiYEnTmTSFxcEiNGXMuCBf2tYtLkKxHhgQeaMn9+P0qUKMhHH/1Ms2YT2b37mNuhGQPz56evjxqV3q88AIgvxo8RkSigLfCtNzrGRkdH6/r16zN/iLO0ycFNiFmxYg9Nm1ZBUv8OBCkR2eCpuA54Wd7DAtzWrX9x003T2LbtCOXKFWHhwgFER1/sdlgm1O3ZA5ekT+/F4cNQunS+h3G+9y+f9PjN1wolS8RMEDt2LI5u3aazceP+tH3NmlUN+kTM+L/atcuyZs0w2ratQaFCEVSpUsLtkIyBqlXhtdfSt5s2dS+W82DlV8b4qdSKyXnzfmHo0E9sFHTjd6KiCrNw4QCWLRtMhQpONVtycgpJSTbMinHRXXdBmTLO+i+/QADcOy0ZM8YPZayYvOKKcnz0UV9rDTN+KSIi7Kwiksce+4oOHT60QWKNe8LD4Ysv0rczrvspS8aM8TOzZm0+q2Jy5cohVKtmFZPG/x0+fIbJkzfx5Zc7PD8mAnM0dBME6tVLX9+92704csmSMWP8yCuvrKJPn9lWMWkCUpkyRVi3bjj161fkt9+O0KjReBYs+NXtsEwoCguDYcPcjiLXAjMZs2EtTJCqWrUkYWHCyy+35Z13OtsckybgVK1akuXLB9O79xWcPJlAly7TeOmlFdbn0bjnkUfcjiBHgZmM2ZyUJohk/EeqV68r2Lr1Lh58sKn1ETMBq2jRAsyY0YtnnolBFR5++EtGjVrpdlgm1Fx0kbM8cgSefNLdWHIQmMlYKhvWwgS4nTuPEh09jtWr96btq1WrjIsRGeMdIsJ//tOSOXP6cPXV5Rky5Bq3QzKh5tFH09efeQbWrnUvlhwEdjJmTABLrZj87rv9PP74ErfDMcYnevSow8aNt1OmTBEAkpJS2Lz5oMtRmZBQuPDZlZSNGkGCf07fZcmYMS7IPMfkRx/1cTskY3wmPDz9n5oHH1zMddeNZerUH12MyISMNm3O7jO2dat7sZyDJWPG5CNV5cUXl9O79yyrmDQhJyVFiY1NJD4+mQEDPuLRR78kJcU69hsfe/rp9PUU/xyQ2JIxY/LRv/61iEce+QrAKiZNyAkLE955pzOvv96B8HDhhRdW0K3bdE6ciHc7NBPMChSAmjWd9cWL3Y0lG5aMGZOP2rSpQbFiBZg9u7dVTJqQJCLcfXdDFi26haioQnz66a80aTKB7duPuB2aCWZxcc7y4YfhkP8NRmzJmDE+Fh+flLbeufNl7Nx5Lz17XuFiRMa4r3XrGqxbN5w6dcqyZcshnnjia7dDMsHs1VfT1597zr04shF4yZgN+GoCyJo1e6lZ83W++WZX2r6yZYu4F5AxfuTSS0uzevUw7rqrAW+/bfd240M9e0J0tLP+2mvuxpKFwEvGbMBXEyBmz95CTMx77N17grffXu92OMb4pRIlCvLGGx0pVcopYklISOaVV1aRkJDscmQm6Iwc6SwvvdTdOLIQeMlYKhvw1fipzBWTw4dfy5Qp3d0Oy5iA8MADi3jggcW0azeFv/4643Y4Jphccomz3L7d3TiyELjJmDF+KDExmeHDP02rmHzppTa8+65VTBqTW7feWo+LLirGN9/spkGDcfz4459uh2SCRalS6ev79rkXRxYsGTPGi/r3/4gJEzZSqFAEs2f35qGHmlnFpDHnoWHDSqxfP4IGDS5m165jNGkygblz/XOgThNgypVLX2/TBgYMgAkT3IsnA0vGjPGiO+64jsqVS/DNN4OsYjIAiMgEEVklIk9k8351EVkgIstEZHRuzzMX5uKLi/PNN4Po3/9qTp9OpHv3GTz77Leo2gCx5gIV8gywvXUrTJ0Kjz/ubjwelowZc4EOH07v19K6dQ22bbuHhg0ruRiRyQ0R6QGEq2oToIaI1MrisBeB/6lqC6CyiMTk8jxzgQoXjuSDD7rzwgutEYH16//AcjFzwbZsgYsvhq5dne0//4TERHdjwpIxYy7I7NlbqFbtNRYt+i1tX6FCES5GZM5DDDDTs74YaJ7FMZcB33nWDwIlc3keIjJCRNaLyPpDfjjIZCAQER5+uDmLF9/KlCndCQuzR/7mAlWv7vQXmzkzfZ8fjMpvyZgxeZCxYvLUqQQWLfK/6hyTo6JAai/eI0CFLI6ZDTwpIl2A9sBXuTwPVR2rqtGqGl0uY18Vc97atKlB8eIFAYiLS6Jbt+msWvW7y1GZgFagQPr61KnuxeFhyZgx5ykxMZkRI86umBw9up3LUZk8OAUU9qwXI4v7oao+C3wGDAPeU9VTuTnP+M7rr69h3rxfiIl5j8mTv3c7HBPIunVzlt98424cBNpNxEbfNy47diyOjh2nMn68VUwGgQ2kP2KsB+zK5rjvgarAK+d5nvGB++5rzD33NCQhIZnBg+fxwAOLSEpKcTssE4j69XOW1au7GweBlozZ6PvGRapKt27T+fLLHZQvX5SlSwdaxWRgmwvcKiKvAH2AzSLybBbHPQS8oqpnsjnPRqDOR5GR4YwZ04GxYzsTGRnGK6+spnPnqRw7Fud2aCbQVKzoLJcvx+3qkMBKxlLZ6PvGBSLCs8/ewLXXXsSaNcNo1Kiy2yGZC6CqJ3A6468GWqnqJlX921AVqvqkqk45x3nH8ydik9Hw4dfx1Ve3UbZsERYt2k6jRuM5dOi022GZQHLllenrce4m84GZjBmTj7ZvP5K23rx5VdatG061aqXOcYYJFKp6VFVnquqB/DjPeFeLFpewbt1w6tatQL16FShbtojbIZlAUqYMFHQKQ6xlzBg/paq89NIKLr/8DebNSx8B3MrrjfEf1aqVYsWKIUye3C2t7+aJE/E2QKzJnfh4Z9mrl6thWDJmTBZSKyYffvhLkpOVnTuPuR2SMSYbxYoVoEiRSABOn04gJmYyQ4d+Qnx8ksuRGb9XrZqz/PFHOHnStTAsGTMmk+PH/14xed99jd0OyxiTC99/f4CtW/9i0qTvueGG9/nzz1Nuh2T82VfOEEXs3QslSsAHH7gShiVjxmSwa9cxmjadaBWTxgSoZs2qsmLFEKpUKcHKlb8THT2O777b73ZYxl+VLg2FC6dvz5/vShiWjBnjkTp0xZYth7jiinJWMWlMgLrmmotYt244TZtWYe/eEzRvPpGZMze7HZbxR6VKwXffwW23OdszZrjSmd+SMWM8RIRx47pw002Xs2LFEKuYNCaAVahQjCVLbmPIkPrExibRt+9sVqzY43ZYxh/Vrg2PPpq+PWVK9sf6iE+SMRGZICKrRORvY/ZkOu4tz1RTkn4AACAASURBVJxvxrhCVVm9em/adoMGlZg372ZKlSrkYlTGGG8oWDCC8eNv4tVXb2TIkPo0bVrF7ZCMv6pdO3190qR8/3ivJ2Mi0gMIV9UmQA0RqZXNcS2Aiqr6qbdjMCY3UismmzSZwIwZP7kdjjHGB0SEe+9tzPjxN6UNfbFz51F27jzqcmTG77z6qrNcuhQO5O8Qgr5oGYsBZnrWF5M+h1saEYkExgG7RKRrVhcRkREisl5E1h86dMgHYZpQlrliMiLCntgbE8xSE7GTJ+O56abpNGw4nm++2eVuUMa/ZBxrbOZMSEzMt4/2xb9ARYF9nvUjQIUsjrkN2AK8BDQUkXsyH6CqY1U1WlWjy5Ur54MwTaiyikljQldKilK5cgn++usMbdpMYezYDW6HZPxFpUpw0UXO+r33wqhR+fbRvkjGTgGpdaLFsvmMa4CxnqlEPgBa+SAOY/5mzZq9NGo03iomjQlRJUsWYv78fjzwQBOSklK4/fb53H33QhITk90OzfiDf/wjfX379nz7WF8kYxtIfzRZD9iVxTG/ATU869HAbh/EYcxZkpJSuO22uRw8eJo2bWpYxaQxISo8PIxRo9oxeXJXChQI580319G+/YccPnzG7dCM2554At55J98/1hfJ2FzgVhF5BegDbBaRZzMdMwFoJSLfAv8A8q8t0ISsiIgwZs3qzT//2ZCFC/tbxaQxIW7gwPosXTqQChWKsmTJTj755Be3QzL+IDzcWU6YkG8fGeHtC6rqCRGJAdoCL3keRW7KdMxJoLe3P9uYzBITk1m4cBtduzply3XrVuC11zq4HJUxxl80aVKF9etHMHXqjwwaVN/tcIw/iIpKXz94EMqX9/lH+qSETFWPqupMTyJmjCtSKya7dZvB5Mnfux2OMcZPVa5cgn//u1laxeWvvx5m9OiVqAsjsRs/0KNH+vqSJfnykVbPb4JS5orJOnXKuh2S8SERqSgiz4vIEyJS3O14TOBKTEyma9fpPPjgF9xyy8fExubf8AbGT3iScgD69cuX6ZEsGTNBZ+3afVYxGXqmAJuBY8BbLsdiAlhkZDgvvNCaokUjmTr1R66/fjL79p1wOyyT32bOzPkYL7JkzASVOXO20LLlZKuYDD0FVPVDVX0DsDlvzAXp2rU2q1YNpVq1Uqxf/wcNGoxj7dp9OZ9ogkfv/O3WbsmYCRrx8Uk88shXxMUlMWzYNVYxGVrKiUh/ERkAlPes9xeR/m4HZgLT1VdXYN264bRseQn795/i+usn8eGHP7gdlglSloyZoFGwYASfftqPV15px9ixXYiMDHc7JJN/ZgC1gJoZ1lO3jcmTsmWLsHjxrdx++3XExycTF5fkdkjGDb/4fsgTrw9tYUx+On48jpkzNzN8+HUA1K5dltq1rbN+CHpBVeNTNzzz396qqhNdjMkEgQIFwnnnnc7ccktdmjevmrZfVdOqL02QGzIEVq706UdYy5gJWKkVkyNGzGfcOJtfLlSJSDjwrYg8LY5BwP1Ad3cjM8EkYyL2008Hadx4Ar/9dsTFiIzP3Xyzszx0yOcfZcmYCUiZKybbtr3U7ZCMS1Q1GYgFtgPdcOa+nQbYMyXjE4899hVr1+6jYcNxfPnlDrfDMb7y3/86y99+g8OHffpRloyZgJOxYrJ16+pWMWkAFNgHLASicKZYy3FwIBGZICKrROSJbN6PEpGFIrJeRN717IsQkT0istTzutqL38MEgA8+6EGXLpdx9Ggc7dt/wJgxa2yA2GBUrlz6etmycPSozz7KkjETMFSVl19eQa9es9IqJj/7bIBVTIY4EemLk3hVAaYD7wIFgEoi0ie7ikoR6QGEq2oToIaI1MrisFuBD1U1GiguItFAXWCaqsZ4Xj/64GsZP1aiREHmzr2Zxx5rTnKycu+9nzNixKckJCS7HZrxprJloWvX9O3Nm332UedMxkQkXERuFJFWGfaJiPTyWUTGZOPMmUQmTnSmNXrxxTZWMWlSVQCqAjVwKihvB4oDhYCLgOxG/I0BUkd2XAw0z+KYw8BVIlIKJ9n7HWgMdBaRtZ6WtSwLoURkhKdFbf2hfOhzYvJXWJjw3HOtmTq1B4UKRTB+/EY6dZpqLWTBZu5cuOYaZ92HUyPl1DI2FegL3CkiY0TkXuAHsr5pGeNTRYsWYMGC/syZ0+eseeRMaFPVMThJ0g7gNDABOA5sV9XXVPWlbE4tivNoE+AITlKX2XLgEuCfwM+e49YBbVS1IRAJdMwmrrGqGq2q0eUyPu4wQaVfv6tZtmwwlSoVZ+DAenZfCkabNjnLLVt89hE5JWNVVHUITkJ2E1AQaKGq9/ksImMy2LXrGM89923ar80aNaLo0aOOy1EZPxQGHAIGAjcCw3JxzimgsGe9GFnfD58E7lDVZ4CtwGDgB1Xd73l/PU5rnAlh0dEXs3Xr3dxyS920ffv3n3QxIuNV997rLGfMgLg4n3xETslYIRFpAjTF+UW4HLhCRJr6JBpjMkitmHziia+ZMGGj2+EYP+V5TFgYaAjsxHn0+BzpiVZ2NpDeyl8P2JXFMVHA1Z7hMxrh9E2bIiL1PPu6AZsu9DuYwFesWIG09Y0b91Or1us8/fRSUlLssWXAa57hYeDnn/vkI3Ia9HUTMDyLdQV8OwKaCWlz5mzhlls+Ji4uidatq9Or1xVuh2T8lKom4SRiqb4XkYeBnjmcOhdYJiIXAx2Am0XkWVXNWFk5EpiE86hyFc6QGatwunAI8Imqfumdb2KCxbp1f3DmTCJPPfUNP/54kPfe60bRogVyPtH4py5d0tfPnPHJR+TUMvYYcACnL8Y/VXWw5zXEJ9GYkJe5YnLoUKuYNDnLPLyEqp5Q1Ukicmd256jqCZxO/KuBVqq6KVMihqquVdUrVbWYqrZV1VOq+pOq1lXVq1X1cV98HxPYRoy4jvnz+1OiREHmzPmZ5s0nsWfPcbfDMnkVGQn9+vn0I3JKxt4HNgPHgLd8GokJeYmJydx++3z+/W+noeGFF1ozbpxVTJpceQNARL7wLBd49t92rpNU9aiqzlTVAz6Oz4SYjh1rsXr1UGrWLM333x8gOnosy5fvcTss46dySsYKqOqHqvoGTlm3MT4TF5fEmjX7KFQoglmzevPww82tMsnkSES6ArEiUhWIzLS0sRSNa+rUKcfatcNo06YGhw6doUuXaZw4EZ/zicZ/jRgBPhi+JKc+Y+U8AyYKUD7j4ImqOtXr0ZiQVrx4QebP78cff5ykUaPshoYyJp2IdMTphB8JvIhT2fgiUNuzrORedMZAVFRhPvtsAA8+uJjmzatSokRBt0MyefHLL87y9GlneqRa3i2izulXYwLOza0mMMOznvoy5oKtXbuPf/3r87ShK6pUKWmJmDkfS3E62ceraj9gk2e53rO0iQON6yIiwnj11fZnFSJ98cV2jh6NdTEqc14+/jh9/fRpr18+p5axI6r6tNc/1RjOrpisX78iAwfWdzskE3jqApOBMiIyEWcYionAdZ7lZSIy0YqOjD9Zu3YfXbpMo2rVknzyST9q1y7rdkgmJ1WrQr166QPAellOLWONReTXTK9tIvKrT6IxISGrisn+/W2uZXP+VHU1UAdn2J1oYA3wAtAO5zFlDDDarfiMyUrFisWoXbss27YdoXHj8Xz22Ta3QzIuyykZW6Oql2V61VLVy/IlOhN0EhOTueMOq5g03qPOM+4I4G6cxOsanNlDflHVrarqu9l9jcmDqlVLsmLFEHr2rMPx4/F07jyN0aNX2ryW/i4pyVlu9P4g5DklY7O9/okmZJ04EU+nTlMZO/Y7q5g0XiMiw4HUfmKrcAaobisiG0TkfnejMyZrRYsWYObM3jz1VEtSUpQHH/yCwYPnEReX5HZoJjubPb/rhgzxekXlOZMxVX3Tq59mQlpERBjHjsVRvnxRvv56oI2qb7xlJTBIVc+IyLXAMVV9GGgDbBcR64xo/FJYmPDkkzHMnt2bIkUi+eijn9m9+5jbYZnsvPpq+npiolcvLYHQLBodHa3r16+H1BaUAIjZZO3AgVPExiZSvXqU26EYPyciG1Q1OpfHLlPVFiIyHxivqnMzvLdEVW/wWaC5kHYPMyYb339/gD//PMWNN9Z0OxRzLql5SHw8FMh+iqvzuX+BDYhofOyjj37mtts+Tpsst2LFYpaIGV+IE5HLgNLAvSKyMHU0fsA3k8kZ40X161c8KxGbNGkj06f/5GJE5px+8u7/G0vGjE+oKqNGraRXr5lMmfIDH330s9shmeBWAngdGAAkAd8AqU3o1pRuAsqvvx7m9tvn06/fHB5//Ku0H7PGj0yc6NXLWTJmvC61YvKhh75A1amY7NmzjtthmSAlIj2AssBMnHHHjAlotWqVZvTodoSHC88/v5zu3Wdw8qRNo+QXbrzRWeZnB35jztfx43F07jwtrWJy5sxeVjFpfK0iznRIxXCSMmtGMAFNRLjnnkZ8/vktREUV4pNPfqFJkwns2HHU7dBMly4+uawlY8Zr9u8/SbNmE1m8eDvlyhXh668H0rv3lW6HZYKcqr4FHAI6Agtw5tK1hMwEvDZtarBmzTBq1y7L5s2HaNBgHMuX73E7LAPw1ltevZwlY8ZroqIKExVVmDp1yrJmzTAaN7Y5Jk2+OQL8A2dqpCJAayC11MkSMxOwatUqw+rVQ+nYsRYpKUr58kXdDim0JSSkr//1l9cum9PclMbkKCVFCQsTChWKYO7cvoSHh1GqVCG3wzKhpYiqbheRxUCYqo7K8F5xt4IyxhtKlizEJ5/czLZtR7jssjKAUySVnKxERFibSr665x643zOW9NKl0KuXVy5r/xdNnqVWTHbvPoPk5BQAypQpYomYcUM7z3ISsDjTe8PzORZjvC48POysCcXfeGMtbdtO4a+/bOSWfBURkT7W2KuvQkqKVy7rk2RMRCaIyCoReSKH4yqIiPcneTI+l7Fi8pNPfmHJkp1uh2RCmKqe9iyPquoPmd77zZ2ojPGNM2cSefnllSxduouGDcfx008H3Q4ptAwZ4ixXrIDh3vmt5/VkzFNmHq6qTYAaIlLrHIePAgp7OwbjW8ePx501x+TMmb1o2/ZSt8MyxpiQUKRIJKtWDSU6+mJ27jxGkyYTmDdvq9thhY6BA9PXvTTemC9axmJwxvsB53FB86wOEpEbgNPAgWzeHyEi60Vk/aFDh3wQpsmL3buP0azZRL74YkfaHJNWMWmMMfmrUqUSfPvtIPr3v5pTpxLo3n0Gzz+/jECY4jDgNW8OY8emb3uhI78vkrGiwD7P+hGgQuYDRKQA8B/gkewuoqpjVTVaVaPLlSvngzDN+frttyM0ajSezZsPUadOWVavHmoVk8YY45LChSP54IPuvPBCawAef3wJ//nP1y5HFQJEzn48mZx8wZf0RTJ2ivRHj8Wy+YxHgLdU1aanDyCXXFKSq64qT+vW1Vm5cqjNMWmMMS4TER5+uDnz5t1MjRpRDB16jdshhY7CnlQn43AXeeSLZGwD6Y8m6wG7sjimDXCXiCwF6ovIeB/EYbxAVYmLSwIgMjKcjz/uy2efDbCKSWOM8SNdulzO1q13pf1IVlV+/fWwy1EFudhYZ/nVVxd8KV8kY3OBW0XkFaAPsFlEns14gKper6oxqhoDfK+qw3wQh7lAqRWTnTtPJTHRaYYtXrwgkZHhLkdmjDEms4z35tGjV3H11W/z3nvfuxhRkAsLO3t5IZe64CtkoqoncDrxrwZaqeomVc12iAtPQmb8TMY5Jpcv38PGjVnWWRgT0HIahkdEokRkoaeY6N3cnmeM2/bsOU5CQjKDBs3jwQcXp40FabxowACvXcon44x5xvqZqar2L3gASq2YTJ1jcunSQTRsWMntsIzxqlwOw3Mr8KGqRgPFRST6PIfvMcYVY8Z04J13OhEREcbo0avo3Hkax4/HuR1WcEmtXP36wosmbAR+c5Z16/alVUzWrm1zTJqgFkPOw/AcBq4SkVJAFeD3XJ5nw/MY191+ezRffnkrZcoU5vPPf6NRo/HWj8ybfvnFWe7bd+7jcsGSMZNm48b9tGw5mT//PM0NN1Rn1SqrmDRBLcdheIDlwCXAP4GfPcfl5jwbnsf4hZYtq7Fu3XCuuqo8v/xymH/8Y4HbIQWPYZ7u7l98AYsWXdClbKJwk6Zu3Qq0bl2D8uWL8PbbnSlQwDrqm6CWm2F4ngTuUNUTInI/MDiX5xnjN6pXj2LlyiH861+LePLJlm6HEzzKlElf79EDTp/O86XsJhLiEhOTOXbM6UcQHh7G7Nm9GT/+JkvETCjIzTA8UcDVIhIONAI0l+cZ41eKFy/I+PE3UaVKScAZ+uLdd9cTH5/kcmQBrH17qFbNWa9b94IuZclYCEutmOzceWraWGIFC0YgqTPSGxPcchyGBxgJjAWOA6WBaVmcZ899TMAZOXI5d9yxgNat3+fPP0+5HU5gKloU3n/fWY+4sAeNloyFqIwVk7/+epidO4+6HZIx+So3w/Co6lpVvVJVi6lqW1U9lcV5x/M7dmMu1I03XkrlyiVYseJ3GjQYx8aN+90OKaRZMhaCsqqYrFPHOhib0JPXYXhs+B4T6K677mLWrRtO48aV+f33EzRvPonZs7e4HVbIsmQsxHz88c9WMWmMMYaKFYuxdOlABg6sx5kzifTuPYsnn/yalBR1O7SQY8lYCFm2bDc9e84kNjaJIUPq2xyTxhgT4goWjGDSpK6MHt2OsDBh8eIdadPfmfxjQ1uEkGbNqtK795XUr1+BRx5pbh31jTHGICLcf38TrrqqPHXrVqBgQUsN8pv9Fw9yJ07EExubSIUKxQgLE6ZP72lJmDHGmL9p1+7StPWUFGXYsE8YPLg+LVpc4mJUocEeUwax1IrJzp2ncfp0AoAlYsYYY3I0efL3TJr0Pa1bv8+4cRvcDsf/LV+ePldlHlgyFqRSKyZ/+ukgp04lcPSoTRBrjDEmd267rR7339+YxMQURoyYzz33LLS+ZFmpnGHu5tjYPF/GkrEglLlicuXKIVSuXMLtsIwxxgSIiIgwRo++kYkTnRlZ3nhjHe3bf8iRI3lPOIJS9epeuYwlY0FEVRk1auXfKiajogrnfLIxxhiTyeDB1/D11wOpUKEoS5bspGHDcezZY+Mcn6Xwhf8ba8lYEFm4cBsPPfQFqvD88zfYHJPGGGMuWNOmVVi3bjjXXFORSpVKULFiMbdDCjpWTRlEOnasxYgR19K6dQ369LnS7XCMMcYEiSpVSrJ8+RBiYxPTfuTHxiZSqJDNZ5zWV+zYMShSJE+XsGQswO3efYywMKFKlZKICO++28XtkIwxxgShIkUiKVIkEoCkpBS6d59BuXJFGTeuC4UKWTrBqlXQs2eeTrXHlAEstWKyU6epnDgR73Y4xhhjQsTmzQdZvnwPH3zwAy1bTuaPP066HZJ7LvWMz3YBLYSWjAWojBWT5coVtbnEjDHG5Jt69SqycuVQLrmkJGvX7qNBg3GsW7fP7bDcUbfuBV/CkrEAo6qMHv33ikmbY9IYY0x+qlu3AuvWDef66y/hjz9O0qLFJD788Ae3wwpIlowFkKSkFO68cwEPPmgVk8YYY9xXrlxRvvjiVkaMuJb4+GRuueVj5s//1e2wAo71uAsgH3/8M+++u4GCBcN5771u9O17ldshGWOMCXEFCoTzzjudqVu3AgsX/kb79jXdDil/paQ4yzNn8nwJaxkLIL16XcEjjzTj668HWiJmjDHGb4gId93VkPnz+xER4aQWf/11hu3bj7gcWT746y9nOXFini9hyZifW7/+j7Q/zCLCyJFtaNKkistRGWOMMX+XOuZYYmIyvXvPokGDcSxZstPlqHysbFlnWTPvLYKWjPmxuXO3cv31k+jUaSpHj9p8YMYYYwJDfHwyxYsX4OjRONq1m8Kbb65FNUir/jt0cJY2tEVwSa2Y7NFjBrGxSTRrVoWiRQu4HZYxxhiTK8WKFeDjj/vyyCPNSE5W7r77M+68cwEJCcluh+aXLBnzM0lJKfzjH1YxaYwxJrCFh4cxcmQbPvywB4UKRfDuuxto23YKhw6ddjs03xg7Ns+nWjLmR06ciKdz56m8845TMTl9ek8efbSFzftljDEmYPXvfzXffjuIiy8uzrff7mbChI1uh+RdyZ7WvkJ5H+/ThrbwI/Pn/8qiRdspW7YIn3xys3XUN8YYExQaNKjEunXDef31NTz0UFO3w/Gujh2dZdGieb6EJWN+pH//q9m//yTdu9ehRo0ot8MxJuiJyATgCmCBqj6bxft3An09m6WANcBdwA7PC+AeVf0xH8I1JqBdfHFxRo5sk7Z94MAppk37kfvuaxzYT4BSY7+AZMweU7ps3rytbN36V9r2Aw80tUTMmHwgIj2AcFVtAtQQkVqZj1HVt1U1RlVjgGXAOKAuMC11vyVixpw/VaV371ncf/9i+vadzZkziW6HdOH27MnzqZaMuSS1YrJ79xl07Pghx4/HuR2SMaEmBpjpWV8MNM/uQBGpBFRQ1fVAY6CziKwVkQkikuUTBhEZISLrRWT9oUOHvBy6MYFNRHjkkWYUL16AWbO20Lz5RH7//bjbYeVNWIZUav/+vF3CS6GcxXODWiUiT2TzfkkR+UxEFovIxyISUuM2ZK6YHD78WkqUKOh2WMaEmqLAPs/6EaDCOY69C3jbs74OaKOqDYFIoGNWJ6jqWFWNVtXocuXKeSlkY4JHp06XsWbNMC69NIqNGw8QHT2OlSt/dzus81e5cvr68bwllF5PxnLT9A8MAF5R1XbAAaC9t+PwV1YxaYzfOAUU9qwXI5v7oYiEAa2ApZ5dP6hq6s/f9UBW9zhjTC7UqVOOtWuH07p1dQ4ePE2rVu8xaVKAVVuKwOWXX9AlfNEyFkMOTf+q+paqfuHZLAcc9EEcfmfPnuM0bz6RRYu2U65cEZtj0hh3bSD9/lQP2JXNcS2ANZo+fPgUEaknIuFAN2CTT6M0JsiVLl2Yzz4bwD33NCQhIZk9ewL0ceUF8EU1Zeam/2uzO1BEmgBRqro6i/dGACMAqlat6oMw89+yZbv58ceD1K5dlgUL+ltHfWPcNRdYJiIXAx2Am0XkWVXN3L3iRuDbDNvPAFMBAT5R1S/zJVpjglhkZDhjxnSgc+fLaNOmhtvh5DtfJGO5bfovDbwO9MzqfVUdC4wFiI6ODooJrQYMqEtiYgpdu15OVFThnE8wxviMqp4QkRigLfCSqh4gi1YuVX0s0/ZPOBWVxhgva9fu0rT1vXtPcMstH/Huu525/PKyLkble754TJlj07+nw/4s4FFV3e2DGPyCqvLaa6vZtOlA2r5Bg+pbImaMn1DVo6o605OIGWP8yOOPL+Gbb3bTqNF4Fi36ze1wfMoXydhc4FYReQXoA2wWkcyDKQ7FeXz5uIgsFZG+mS8S6FIrJu+7bxFdukwLjjFUjDHGmHzy5psd6d69NsePx9Ox41ReeWUV6V03g4vXkzFVPYHTiX810EpVN2Xug+EZSDEqw6CJM7wdh5syV0y+/HJbihSJdDssY4wxJmAUK1aA2bP78OSTLUlJUR54YDFDhnxCfHyS26F5nU/GGQvlpv+MFZNlyxZhyRKrmDTGGGPyIixMeOqpGGbN6k2RIpFMnvw9rVu/T1JSituheZXNTelF69f/QZcu0zhw4JRVTBpjjDFe0qvXFdSsWZquXafTpctlREQE1wRClox50Y4dRzlw4BStWlVjzpw+1lHfGGOM8ZL69SuyadMdlCyZPmPNkSOxlC4d+P/WWjLmRX36XEmRIpG0a3cpBQqEux2OMcYYE1RKlSqUtr5r1zEaNhzHiBHX8cwzrQgLC9yZbIKrnS+fJSWlcN99n7N27b60fZ07X2aJmDHGGONjK1bs4fDhWJ57bhk9eszg5Ml4t0PKM0vG8ujEiXi6dJnGa6+toU+fWSQkJLsdkjHGGBMyBgyoy2efDaBUqULMm/cLTZtOZOfOo26HlSeWjOVBasXk55//RtmyRZg6tae1hhljjDH5rF27S1mzZhiXX16Gn346SIMG4/jmm11uh3XeLBk7T+vX/0GjRuP58ceDXH55GVavHkrTplXcDssYY4wJSZddVobVq4fRvn1NDh+OpX37D9m//6TbYZ0X68B/HubN20q/fnOIjU2yikljjDHGT5QqVYj58/vxyCNfUrlyCS66qLjbIZ0XS8bOU1xcEoMG1efddzvbo0ljjDHGT4SHh/Hyy+3O2rdmzV5q1ixNmTJFXIoqdywZOw9du9Zm9ephNGhwMSKBW0JrjDHGBLtt2w7TocOHlCpViE8/7ceVV5Z3O6RsWZ+xczhxIp5u3aazbNnutH0NG1ayRMwYY4zxc4ULR1KjRhQ7dx6jceMJfPrpL26HlC1LxrKRWjE5b94vjBgxn+Tk4JoHyxhjjAlmlSuX4NtvB3PzzVdx6lQCXbtOZ+TIZaiq26H9jSVjWdiw4eyKyfnz+xEebv+pjDHGmEBSpEgkU6f24LnnbkAVHntsCQMGfERsbKLboZ3FMoxM5s3byvXXT06bY3LVqqFcemlpt8MyxhhjTB6ICI891oK5c/tSrFgBZs3awg8//Ol2WGexDvwZvPXWOu6+eyGqWMWkMcYYE0S6dq3NypVD+PHHgzRqVNntcM5iLWMZ1KxZmvDwMJ577gYmTrzJEjFjjDEmiFx9dQX69786bXvBgl+ZMmXThV/4+HFnuX17nk4P+ZaxlBRNm+m9XbtL+fXXu6lePcrlqHwnMTGRvXv3EhcX53YoxgBQqFAhKleuTGRkpNuh5Cv7uxjcQvXPdSDZv/8k/frN4eTJBH788SAjR7bOe//wAwec5aZN0KnTeZ8e0snYnj3H6d59Bi+80Jq2bS8FCOpEDGDv3r0UL16catWq2RAdxnWqyuHDh9m7dy/Vq1d3O5x8ZX8Xg1co/7kOJBddVJwXX2zDP//5OS+/rvL9tgAAGRpJREFUvJKffjrItGk9KVmy0PlfrEcP+OgjKFo0T7GE7GPK1Dkmv/tuP0888bVflrr6QlxcHGXKlLGbv/ELIkKZMmVCsnXI/i4Gr1D+cx1o7ryzAYsX30Lp0oX57LPfaNx4Atu2HT7/C1W+sD5oIZmMzZ27leuvn8SBA6eIianG558PCKkbYih9V+P/QvnPYyh/92Bn/28DR6tW1Vm3bjhXXVWerVv/omHD8SxZsjNfYwipZExVeeWVVfToMYPYWGeOyUWLbrHJvs0FSUhIYN++fW6H4bd27NjhdgjGJfPmzWPt2rVuh+ET9uc6uNSoEcXKlUO46abLiY1NpEiR/O3rF1LJ2GOPfcUDDyxGFauYNF7z0EMP8fvvv7sdhiv279/Pl19+ycmTJ7M95u233+bbb7/Nx6hyT0QmiMgqEXkim/fvFJGlntf3IvJubs4LJJ9++imrVq1K2+7fvz/x8fFnHZOUlJTt+anv/fHHH9x2220MGjSIW2+9lV27dnH33XczduxYHnzwwbP+jiQmpg+4+dprrzF9+vRcfZY/+fTTT5k6darbYRgvKl68IB9/3JeVK4fSuHH6Y8eUFN93YwqpZKxDh1qUKFGQadN68thjLawZ2QWDBg3immuuoUmTJvTu3ZvExEQOHDhA27ZtadKkCf/+978BSE5OZsSIEbRo0YKBAweSkuKf01GtWrWK8uXL07hxY7dDOUtiYiJdunShWbNmTJw4MdvjvvvuO9q0aUOzZs0YPXo04Pzib926NfXr1+d///sfAHv27CEmJoYbbriBESNGoKr8+uuv9O3blxUrVtCyZUsSEhKy/NyRI0cyZswYv/t/KCI9gHBVbQLUEJFamY9R1bdVNUZVY4BlwLjcnBdImjRpwksvvQTAihUriIqKomDBgmcdc9ddd9GyZUtiYmKoX78+VatWJSYmhpYtWzJ48GAALrroIsaMGcPRo0eZPHkyDz74IG+++SZjxoxhy5YtVM7Qp6Zp06Zcf/31xMTE8Nprr/HOO++kXa9JkybExsZe8PcaOnQoTZo04dlnn832mJ07d9KpUydatGjBAw88AMDx48fp0KED7dq1o3v37iQkJGS5795772XBggXn/CFiAk9YmHDttRelbc+Zs4VWrd77//buPizKKm/g+Pen+FJLiWllL75UWm66Wqllaga6udSzZl6Zmpr0orllmPlUZpqZadKbtdolD6ahtqnY2tqaqYkyCYqt2LPa5dv2rGBhZaSIC5gK/J4/bhwFZmRQYIbh97muuWRmzj3n95u553juc5+5Dz//nFe19VbpqweAM5c86NmzJRkZTzN4cHs/RmRmz55NamoqoaGhJCYm8vLLL9O/f382b95MSkoKmzdvJiEhgePHj5OcnEyzZs1YsWKFv8P2aNGiRTz11FPVVl+/fv0IDw933+bOneux3OzZs+nUqRObNm3ir3/9q9f/MKKjo4mPjyclJYXly5eTnp7Oe++9x9SpU/nnP//J2rVrycrKIi4ujtjYWDZs2MD333/PN998w44dO4iPj+fll1/m2muvJT093WO9ISEh9O7dm02bNlXlW3MuwoFlxX9/AfTwVlBErgIuV9U0X7cTkcdFJE1E0rKysior5kqVmZlJ+/btyc7OJjw8nIkTJ7Jz507CwsJKjFDFxcXx5Zdf4nK5ePfddxkyZAgul4svv/ySDz/8EHDmSKWmpnLXXXdRUFBAly5dmDlzJrfffju7d+8mIiLCfcoyMjKS5557jqioKIYMGcLw4cOZOHEijzzyCJGRkVxwgeepI6NGjSqx/0+dOtVjuU8++YTCwkJSU1PZt28f3377rcdy48eP56WXXiI5OZnMzExcLhcfffQR48aN44svvqBZs2asWbPG42MAw4YNC9i2yZy/goIiJk7cwMaN++nS5X22b/+pyuoK6ktbbNv2A/fdl8CcOffQt+8NADY/7ExVNTLowy9TVZXc3Fzq16/PV199RXR0NCJCt27d2Lp1K19//TX/VXytlkGDBpGX5/mo5JdffuHhhx/m0KFDXH/99cTHx/Poo48yZcoUWrVqxZQpUwgPD6dVq1ZMnDiR+vXrAxAfH8/IkSMZPXo0N910E6NGjWLEiBG0aNGChx9+mJycHPr27cuECRPOmsexY8do1KgR4JymGThwICJCz549mT59OhkZGWXqPXjwYJk6PG3ryaefflruewvgcrmIiYkBoGfPnqSlpREREVGm3OHDh2nevDkATZo04ejRozRp0oQdO3bQunVrjh8/TlhYWIl4Dh06RNOmTenQoQMFBQWsWrWK7OxsWrdu7bXerl27kpyczB133OFT/NXkN8CpyX6HgVvOUnY0EFuR7VR1LjAXoHPnzmf/Uvjpu1i/fn0iIyNZsGBBicfDw8OpW7fkFI4xY8Ywa9asMq9RWFjoLjt//nxUlaSkJN544w2OHTtGs2bNaNu2LWlpaeTm5gLw4osvkpGRwauvvkpeXh6tW7cmISGBzz//nClTpniNNy4uzoeknf1/4MCBAPTp04eUlBTatCk7gPmvf/2LW25xPr7LLruMnJwcnnzySffzWVlZXHbZZdx7771lHgPo2rUrr776Kg899JBPcZmaJSSkDklJUfTvn8BXXx2gW7cPWLToPu6//8ZKrytoR8ZOrTGZmXmUuXO/rjWXrqgJoqOjadWqFZdffjm9evXiP//5D78pvjbLhRdeyNGjRzl48CCXXOKsCXrLLbd4/U/8tddeY+jQoaSmptKuXTv279/vtd6VK1cyatQo4uPjARgwYACrV68GYM+ePXTp0oUZM2YwaNAgNm/ezIoVKzh0yPefOB84cICYmBhWr17NypUrvdbrqQ5v256rvLw8rrrqKgAuueQSDh70vA5b9+7dee+991i8eDEZGRl06NCByMhItmzZwqxZs+jVqxchIaeP2RISEmjXrh1XXnklALm5uSxbtoyWLVsiIl7rveCCCyrl1FMlywVOHZ2F4qU9FJE6QATgqsh2NUFRURFpaWklRpvCw8P56aefKCwsLFHW5XK5/168eLH7tOJrr70GQGJiIn//+995/fXXEREKCwtZsGABM2fOZOzYscTGxlKnjvNW7dmzh4ULF/L6669z88038+yzzzJy5EiSkpJYuHDheefl6/4/YMAAXnnlFVauXMmaNWvo3bu3+7nU1FSys7NLTEEo/ViA7temEl1xxUW4XA8zfHhH8vNPMmDAx7zyiqvS55EF3ciYqvLnP3/FuHFrS6wxafPDPPBTB3X27NmkpKTQoEEDRISLL77YfcScl5dHixYtSjy2YsUKcnNzGTZsWJnX2rNnD6NHjwacifSlndlQ9unTp0TD2rt3b2JjY9m9ezedO3cGYO/evaSmprJgwQLy8vL44YcfaNKkiddcioqKUFVEhJCQEF555RVCQ0NLnBYsXa+nOrxtW1q/fv3IObXsBs5k68cff7xMudDQUPeoXW5uLqGhoR5fLy4ujqSkJCZPnsz48eMREWJiYli2bBkiwpgxY1i3bh19+vRh3759vPXWWyQmJrq3DwsLY+HChTz00ENs3brVa73p6enuEbgAsg3nFOMWoCOw10u5O4Cv9PQRna/b+c5P38UGDRqQlpZGXl4ey5cvZ+DAgYSFhbF//353x+mUM+8PGTLEPQJ6SosWLdzzxwDy8/OJjIykadOmtGrVil27dgHOnKy4uDimT5/OiBEjyMnJ4a233qJevXpMmzaNJUuWcOLECfdo8plGjRrF3r2n3+5evXoxefLkMuVO7YfgHDB4m684adIkUlJSePPNN4mKinLvr4cPHyY6Oprly5e7y3p6LED3a1PJGjYMYcGCfnTocBnPP5/IlClf8ssv+cyefU+l1VFjj+g8KSgo4qmnPueZZ5yO2LRpEfaLyQA1atQo5s+fT2FhIbfddhsulwtVZdOmTdx66610796ddevWAbBu3TrCwsI8vk7btm3ZunUrAI8//jiJiYnUr1+frKwsCgsL3a8BlOmQhISEcOmll7JkyRIGDBgAwA033EBMTAwul4sXXnjBPTrnTZcuXdy/FJw5cyYTJkxg3rx5JTr/pev1VIe3bUv79NNPcblc7punjhhAp06dSElJAWD79u20atXKY7m6detyww3OKfyhQ4cCzn8w33//Pb/++itff/01IkJ2djYPPvggH3zwgfu07BNPPOHO/ciRI4SFhXmtd8WKFfTp08drXn6yAnhIRGYCA4GdIuJptvcfgI1n2W5VlUdaRQYPHkxKSgozZ87kuuuuY+zYsQCsXbvW677lzfXXX0/jxo0ZOnQoaWlppKenc+GFF5Kfn09ubq67c3RqbuJjjz2GqpKVlcUPP/xAQUEBTz/9NMnJySxatMhjHXFxcSX2f08dMfB9/we46aab+O677xg3bhzgXKrmgQceYMaMGbRs2dLrYwDLli3jj3/8Y4XeJ1MziQj//d/d+OyzB7niilAeeeTmyq1AVQP+1qlTJ1XnsNS5eREV9TeFKdqgwau6ZMk3XsvVZrt27fJr/VFRUZqcnKyqqtHR0ZqQkKA//vij9u7dW2+77TZ99tlnVVU1Pz9fBw8erN27d9dhw4ZpYWGhx9f7+eef9e6779aePXvqiBEjtKioSNevX68RERE6cuRIfeCBBzQpKUnT09M1KiqqzParVq3SNm3aaFFRkaqq/vjjj3rPPfdot27ddOjQoXry5Mmz5pOfn6/9+/fXvLw8Xbp0qbZv314jIiK0Xbt2mpmZ6bFeT3V42vZ8ZGRk6I033qhjxozRzp07a0FBga5fv15nz55dpuzw4cN148aN7vufffaZXnPNNRoaGqqDBw/WgoICff7557VZs2Z655136p133qkul0v37dun3bt31x49eujUqVO91pucnKwTJ048a7ye9ksgTau4bQEa43SomlXldu42rJycq9OBAwe0b9++unXrVp08ebKqqubk5Gh6errOmzdPJ02apNu2bdOTJ09qYWGhduzYUVVVk5KSdPz48e7XOX78uPv789xzz+m3336r48aN06ioKN21a5fGxsZqUlKSvvHGG5qUlFQmjnfeeUeXLFlSqbnl5ORohw4d9JlnntG2bdvqkSNHdOfOnR73w8mTJ+uiRYvc9+fMmaNhYWHufX3p0qUeH9u7d6+OHDnyrHH4+zM2VSM//0SJ+/v3H1EdM8bpn7z7rqpWvP3ye0fLl5uvnbEtW77X5s1nakrKfp/e0NrIGofKt337dl2zZo2/wyjjwIEDmpCQoEeOHPFrvW+//baeOHHirNv4qzNWXbdA7IxlZmbqpk2bNCMjQ5s3b+7uaLRp06ZE52TZsmXao0cP9/Olb127dtV///vfqqqamJio2dnZevDgQZ0xY4aqqu7YsUN3796tgwYN0v37y7bNMTEx+uGHH1Z6focPH3Yf7FWF2NhYzcnJOWsZf3/Gpup99NEOrV//VZ3X+9nz6oyJs01g69y5s6alpZ3+xdEZMWdl5XHppacX5jx+vIAGDYJuKlyl2b17N7/97W/9HYYxJXjaL0Vkm6p29lNIlcrdhp3BvovBzz7j4DdhQiIxMc5le8awhbdHtCDk/bgKt181es7Yp5/u4dprZ/Hxxzvdj1lHrHw1oQNuao/avD/W5tyDnX22tcOMGb9n/vx7qVdHWczv+Gne0vI38qBGdsZUlXfeSaV//wRyc0+wfn31LuhZkzVs2JBDhw5ZQ2ECgqpy6NAhGjZs6O9Qqp19F4NXbd6va6NHH72ZpLd/x/JrtnL1my/BOaw2UuOGkQqow9NPfc6cOc6Q/7RpEbz4YkBdSDKgXX311WRmZhKoVwQ3tU/Dhg1LLJVTW9h3MbjV1v26tuo+9n4Ye/85b18lnTERmQ/cCKxSVY8Lg/lSprSjNGAwA1g9J40GDeqyYMF9trRRBdWrV49rrrnG32EYU+vZd9EYc0qln6b0ZRHdc11o90HuZzVtaNr0QjZsiLKOmDHGGGNqvKqYMxZO+YvollvG0yK709lAZw6wZctjdOtmVz02xhhjTM1XFacpfVlEt9wy6mGR3Zv0R/6haksbGWOMMSZoVEVnzJdFdCu00O62bdt+EZFTK0A3BX6phDj9KRhyAMsj0ARDHmfm0PJsBWuSM9qwYPiMwPIIJMGQAwRfHhVqv6qiM+bLIroVWmhXVS899beIpNX0C0EGQw5geQSaYMgjGHLw5FQbFiz5WR6BIxhyAMujKjpjK4BkEbkSuBsYLCLTVHXSWcp0rYI4jDHGGGMCXqVP4FfVozgT9LcAEaq6vVRHzFOZnMqOwxhjjDGmJqiS64ypajanfy15zmW8mHtOQQWWYMgBLI9AEwx5BEMOZxMs+VkegSMYcoBankeNWCjcGGOMMSZY1ci1KY0xxhhjgoV1xowxxhhj/Mg6Y8YYY4wxfhSwnTERmS8iqSIy6XzK+FN58YlIIxFZLSJfiMjfRKR+dcfoC1/fZxG5XET+t7riqqgK5DFHRPpWV1wV4cM+1VhEPi9eSiyuuuOriOL9Jfksz9cTkZUisklEHq3O2M5XMLRfEBxtmLVfgSVY2rDKbr8CsjNWlYuNVxcf4xsKzFTVPsBPQGR1xuiLCr7Pb3F6ZYWA4mseInIH0ExVV1ZrgD7wMYeHgI+KLzp4kYgE5EUURaQxsBBnaTRvooFtqtodGCAiF1VLcOcpGNovCI42zNqvwBIsbVhVtF8B2RmjkhYb97NwyolPVeeo6rriu5cCP1dPaBUSjg/vs4j0AvJwGuRAFE75i9PXA94HMkSkX/WF5rNwyv8sDgHtRSQMaA58Xz2hVVghMAg4epYy4ZzOdyMQcI2yF+HU/PYLgqMNC8far0ASTnC0YZXefgVqZ6z0QuKXn2MZf/I5PhG5HWisqluqI7AKKjeP4lMTLwEvVGNcFeXL5zEc2AW8AdwqItHVFJuvfMkhBWdNtDHA7uJyAUdVj/pwsedA/457EwztFwRHG2btV2AJijasKtqvQO2MVfpi437gU3wicgkwGwjUOTG+5PECMEdVj1RbVBXnSx43A3NV9SfgL0BENcXmK19yeBn4k6pOBfYAj1RTbFUh0L/j3gRD+wXB0YZZ+xVYalMbVqHveCA2AHB6IXFwFhLPOMcy/lRufMVHZB8DE1R1f/WFViG+vM+/B0aLiAu4SUTmVU9oFeJLHv8HXFv8d2cg0D4TX3JoDPxOROoCtwE1+arOgf4d9yYY2i8IjjbM2q/AUpvasIp9x1U14G7AxcB2YCbOMGVHYFo5ZRr5O+5zyOEJIBtwFd8G+Tvuc8mjVHmXv2M+j8/jIpz/WDYCqcBV/o77HHK4FdiJc1S2Dgj1d9zl5OQq/rcX8FSp51oW5/JnYCvOxF+/x1xJn1NAt18VyCOg2zBrv/wf+znkUWPasMpsvwJ2OaTiXyvcBWxUZ9j1nMr4U6DH5yvLI3AEQw4VISJX4hxdrtXy52gEjGBov6BmxFieYMgBLI+aqCLtV8B2xowxxhhjaoNAnTNmjDHGGFMrWGfMGGOMMcaPQvwdgDEiMgXnAnoHix9qD2QV307iXJH5NZzJnr8CmcAQVT1Z7cEaY4wXInIFcERVj/k7FlOz2MiYCRTTVTVcVcOB94rv9wTicZaVAIhWZxmNXJyfoxtjTLUQkdEicv9Znq8L/A34wMvz9c74+2kRGXzGfRsYqeVsBzCBrjHgPsoUEcG5gN4Jv0VkjAl6IhKKs2TPMZx26FKgQETG4yyDcwHwB1XNFZE6QBywwNlUpqrq5FIvuVlEjgFFQAvgOxH5EyDAhSLSE4jFzgDUStYZM4FiooiMwFnO4+fi+4/jrEv2JM61WmYDlwArgQ3+CtQYE/xUNRfoJiKtcTpJH+KMyg/EGaXfAyAizYuf36Sq/1P82AQRWQfE4FyLqhBYA/wDaApcB+zDad+uAq5T1WPOsSbRqpoiIvE4ZwBWV1fOxn+sM2YCxXRV/Qu455C57xc/Bs7pyh7AcbVrshhjqlDxMk+PAH2AKOAe4BdgFPCRiGQDM4D5OKNYkSLyhzNeIhR4Fueq+O/jzHtthbMO5m9wrpg/qPh1p5Sq284A1DLWGTM1TRyQLCKzio82jTGmKjTHmVd9EbCo1HMngHk4o1sdcQ4Qi0QkEuiqqlNEpK6qFhZ3rADa4nS+xgOP4bRl7wNHgDBOzzWzMwC1kHXGTI2iqtkisgG4H1jm73iMMcFJVbcD20XkLsp2isJV9RMAEWkGLBaRIqAJ0EhEejhPycvFpxwb4YyoTcTpxDXCGTU7CUwCHixe5xPsDECtZFfgN8YYY7wQkY2UPV1YR1V7FT8vOOsOFpQaGauD839soYg8ijMqdurHSNfhTNI/UHy/Hs4vMXvgdNZ2AslARzsDUDvYyJgxxhjjRfElds72vAIFHh4vOuPvDzjjkhciMhb4SVWXnrlN8YianQGohWxkzBhjjClFRNrgzOnK9fQ0ziT8cTiT+AEUZ67XxUAGzmhXiKre6eG1xwMHzvyRkqndrDNmjDHGVAERaaCqx/0dhwl81hkzxhhjjPEjWw7JGGOMMcaPrDNmjDHGGONH1hkzxhhjjPEj64wZY4wxxvjR/wNZqCjwrwgqwwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "X=data.loc[:,['PM2.5','PM10','SO2','CO','NO2','O3']]\n",
    "Y=data.loc[:,'有无污染']\n",
    "modelLR=LM.LogisticRegression()\n",
    "modelLR.fit(X,Y)\n",
    "print('训练误差:',1-modelLR.score(X,Y))  #print(accuracy_score(Y,modelLR.predict(X)))\n",
    "print('混淆矩阵:\\n',confusion_matrix(Y,modelLR.predict(X)))\n",
    "print('F1-score:',f1_score(Y,modelLR.predict(X),pos_label=1))\n",
    "fpr,tpr,thresholds = roc_curve(Y,modelLR.predict_proba(X)[:,1],pos_label=1) ###计算fpr和tpr\n",
    "roc_auc = auc(fpr,tpr) ###计算auc的值\n",
    "print('AUC:',roc_auc)\n",
    "print('总正确率',accuracy_score(Y,modelLR.predict(X)))\n",
    "fig,axes=plt.subplots(nrows=1,ncols=2,figsize=(10,4))\n",
    "axes[0].plot(fpr, tpr, color='r',linewidth=2, label='ROC curve (area = %0.5f)' % roc_auc) \n",
    "axes[0].plot([0, 1], [0, 1], color='navy', linewidth=2, linestyle='--')\n",
    "axes[0].set_xlim([-0.01, 1.01])\n",
    "axes[0].set_ylim([-0.01, 1.01])\n",
    "axes[0].set_xlabel('FPR')\n",
    "axes[0].set_ylabel('TPR')\n",
    "axes[0].set_title('ROC曲线')\n",
    "axes[0].legend(loc=\"lower right\")\n",
    "\n",
    "pre, rec, thresholds = precision_recall_curve(Y,modelLR.predict_proba(X)[:,1],pos_label=1)\n",
    "axes[1].plot(rec, pre, color='r',linewidth=2, label='总正确率 = %0.3f)' % accuracy_score(Y,modelLR.predict(X))) \n",
    "axes[1].plot([0,1],[1,pre.min()],color='navy', linewidth=2, linestyle='--')\n",
    "axes[1].set_xlim([-0.01, 1.01])\n",
    "axes[1].set_ylim([pre.min()-0.01, 1.01])\n",
    "axes[1].set_xlabel('查全率R')\n",
    "axes[1].set_ylabel('查准率P')\n",
    "axes[1].set_title('P-R曲线')\n",
    "axes[1].legend(loc='lower left')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "说明： 1、这里利用空气质量监测数据，建立Logistic回归模型对是否有污染进行分类预测。其中的输入变量包括PM2.5,PM10,SO2,CO,NO2,O3污染物浓度，是否有污染为二分类的输出变量（1为有污染，0为无污染）。进一步，对模型进行评价，涉及ROC曲线、AUC值以及F1分数等。需引用sklearn.metrics中的confusion_matrix,f1_score,roc_curve, auc，以及import accuracy_score等。 2、modelLR.score(X,Y)为预测模型的精度得分（基于训练集的）。分类预测的精度得分为总的预测正确率。也可通过accuracy_score函数得到同样结果。 3、confusion_matrix(Y,modelLR.predict(X))：计算模型的混淆矩阵。 4、f1_score(Y,modelLR.predict(X),pos_label=1)：针对1类计算F1得分。 5、modelLR.predict_proba(X)中存储模型预测为0类和1类的概率，这里关心预测为1类的概率。 6、roc_curve：计算预测为1类的概率从大到小过程中的TPR和FPR。auc(fpr,tpr) 计算ROC曲线下的面积。 7、precision_recall_curve：计算预测为1类的概率从大到小过程中的查准率P和查全率R. 8、ROC曲线和AUC值，以及P-R曲线均表明，该预测模型的预测误差（训练误差）很小。"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
